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In  recent  years,  vegetated  floodplains  and  wet¬ 
lands  have  been  regarded  as  constituents  of  the 
ecosystem  where  significant  transport  processes 
take  place  during  floods.  In  this  sense,  sedimen¬ 
tation  has  been  identified  as  a  major  contributor  to 
nonpoint  source  pollution.  Engineering  tools  are 
thus  needed  for  estimating  both  the  mean  flow  and 
turbulence  structure  as  well  as  the  suspended  sedi¬ 
ment  transport  capacity  of  vegetated  waterways. 

RESEARCH  OBJECTIVE: 

The  two-equation  turbulence  model  based  on  the 
k-e  closure  scheme  was  developed  to  simulate  the 
flow  and  turbulence  characteristics  of  open-chan¬ 
nel  flows  through  nonemergent  vegetation.  Once 
the  performance  of  the  model  was  verified,  the 
flow  structure  of  vegetated  open  channels  was 
num^cally  simulated.  Simulated  rigid  and  flex¬ 
ible  plants  were  used  to  validate  die  model.  Fi¬ 
nally,  dimensional  analysis  allowed  identification 
of  the  dimensionless  parameters  that  govern  sus¬ 
pended  sediment  transport  processes  in  the  pres¬ 
ence  of  vegetation,  and  thus  helped  in  the  design 
of  numerical  experiments  to  investigate  the  role  of 
different  flow  properties,  sediment  characteristics, 
and  vegetation  parameters  upon  the  transport 
capacity. 

SUMMARY: 

The  two-equation  turbulence  model  was  found  to 
accurately  represent  the  mean  flow  and  turbulence 


structure  of  open  channels  through  simulated 
vegetation,  thus  providing  the  necessary  informa¬ 
tion  to  estimate  suspended  sediment  transport 
processes.  A  reduction  of  the  averaged  stream- 
wise  momenrnm  transfer  toward  the  bed  (i.e., 
shear  stress)  induced  by  the  vegetation  was  iden¬ 
tified  as  the  main  reason  for  lower  suspended 
sediment  transport  capacities  in  vegetated  water¬ 
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were  used  to  solve  the  sediment  diffusion  equa¬ 
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concentration  slightly  in  excess  of  the  ones  pre¬ 
dicted  by  the  Rousean  formula.  A  power  law  was 
found  to  provide  a  very  good  collapse  of  all  the 
numerically  generated  data  for  suspended  sedi¬ 
ment  transport  rates  in  vegetated  chaimels. 
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1  Introduction 


Background 

Historically,  vegetation  in  streams  and  rivers  has  been  considered  by 
hydraulic  engineers  as  a  source  of  flow  resistance,  and  as  such  it  has  usually 
been  eliminated  with  the  goal  of  improving  water  conveyance.  This  explains 
why  earlier  research  interests  were  focused  primarily  on  the  estimation  of 
resistance  laws,  mean  velocity  distributions,  and  the  determination  of 
approximate  rules  for  the  partition  of  the  total  action  of  gravity  between  friction 
drag  due  to  bed  roughness  and  form  drag  due  to  plants.  In  recent  years,  however, 
plants  in  aquatic  environments  have  reached  a  different  status,  and  vegetation  is 
no  longer  regarded  merely  as  an  obstruction  to  the  movement  of  water,  but  rather 
as  a  means  of  providing  stabilization  of  banks  and  channels  and  habitat  and  food 
for  animals,  as  well  as  pleasing  landscapes  for  recreational  use  (Haslam  and 
Wolseley  1981).  The  preservation  of  vegetation  is  nowadays  considered  of  great 
relevance  for  the  ecology  of  rivers.  In  recent  years,  such  unprecedented 
environmental  concerns  have  motivated  the  onset  of  several  studies 
concentrating  on  the  characterization  of  turbulent  transport  processes  in  natural 
flow  conditions.  In  particular,  there  has  been  an  increasing  need  for  the 
understanding  of  retention  processes  in  wetlands,  by  which  suspended  solids 
and/or  chemical  contaminants  (pesticides,  heavy  metals,  etc.)  are  being  deposited 
and  retained  within  a  natural  or  artificial  waterway.  In  particular,  vegetated 
floodplains  and  wetlands  have  been  regarded  lately  as  constituents  of  the 
ecosystem  where  significant  transport  processes  take  place  during  floods.  In  this 
sense,  sedimentation  has  been  identified  as  a  major  source  of  nonpoint  source 
pollution  impairment  in  U.S.  rivers  and  Isikes,  where  excessive  sedimentation 
results  in  the  destruction  of  fish  habitat,  decreased  recreational  use,  and  loss  of 
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water  storage  capacity  (U.S.  Environmental  Protection  Agency  1993).  In  Illinois 
for  instance,  field  studies  conducted  in  riverine  wetlands  have  indicated 
sediment-trapping  efficiencies  ranging  from  60  to  85  percent  (Demissie  1990), 
and  estimates  by  the  U.S.  Department  of  Agriculture  (USD A)  indicate  that 
annual  offside  costs  of  sediment  derived  from  cropland  erosion  alone  are  of  the 
order  of  $2  to  $6  billion,  with  an  additional  $1  billion  arising  from  loss  in 
compared  productivity  (USDA  1987). 

All  these  have  prompted  the  development  of  engineering  tools  for  the 
estimation  of  contaminant  and  sediment  transport,  for  the  assessment  of 
environmental  impacts,  for  the  evaluation  of  design  alternatives,  and  for  the 
management  of  wetlands.  Although  some  general  models  based  on  rather  crude 
assumptions  have  been  developed,  they  have  benefited  veiy  little  from  advances 
in  the  knowledge  of  turbulence  in  the  presence  of  vegetation  from  other  research 
areas  such  as  atmospheric  sciences.  As  a  result,  very  few  physically  based 
models  exist  to  help  engineers  evaluate  transport  processes  and  in  particular  the 
sediment  retention  capabilities  of  vegetated  waterways. 


Purpose  of  Study 

The  overall  objective  of  the  present  work  is  to  investigate  the  effect  of 
vegetation  on  the  mean  flow  properties  and  on  the  turbulence  stmcture  in 
open-channel  flows  and  the  implications  of  the  resulting  flow  structure  for  the 
entrainment,  transport,  and  deposition  of  suspended  sediment.  The  working 
hypothesis  is  that  if  vegetation-induced  roughness  increases  flow  resistance  via 
momentum  diffusion,  the  same  roughness  should  also  reduce  the  diffusion  of 
suspended  sediment.  In  particular,  it  is  important  to  know  if  the  suspended 
sediment  distribution  is  significantly  different  from  the  Rousean  distribution  for 
equilibrium  open-channel  suspensions,  and  if  it  is,  what  are  the  implications  for 
the  advective  and  diffusive  transport  of  sediment  in  vegetated  channels. 


To  achieve  such  objectives,  knowledge  about  turbulence  characteristics  in  the 
presence  of  vegetation  coming  from  atmospheric  boundary  layers  will  be 
coupled  with  advances  in  the  numerical  simulation  of  free-surface  flows.  This 
integration  will  be  used  to  produce  a  two-equation  turbulence  closure  scheme  for 
modeling  the  complex  turbulence  stracture  of  flows  through  vegetation  and  for 
estimating  related  transport  processes  in  plant  environments.  The  numerical 
model  will  then  provide  quantitative  information  about  the  role  played  by  flow 
parameters,  sediment  properties,  and  vegetation  characteristics  in  the  suspended 
sediment  transport  capacity  of  vegetated  free-surface  flows,  and  therefore  will 
facilitate  the  assessment  of  the  factors  that  influence  the  sediment-trapping 
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ability  of  wetlands,  as  well  as  the  conditions  under  which  previously  deposited 
sediment/pollutants  might  be  reentrained  into  suspension  and  exported  out  of  a 
given  system. 
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2  Literature  Review 


Mean  flow  and  turbulence  characteristics  in  the  presence  of  vegetation  have 
received  a  lot  of  attention  in  the  last  few  years,  specially  for  the  case  of 
atmospheric  flows  over  plant  canopies.  One  of  the  main  motivations  for  such 
studies  has  been  the  need  for  understanding  related  transport  processes  in  natural 
environments,  such  as  the  transport  of  pollutants,  heat,  carbon  dioxide,  etc. 
Regarding  free-surface  flows  in  streams,  engineering  research  on  vegetated 
open-channel  flows  has  traditionally  been  limited  to  the  estimation  of  resistance 
laws.  In  general,  the  investigations  may  be  classified  into  two  groups 
corresponding  to  the  study  of  rigid  and  flexible  vegetation,  respectively.  An 
extensive  bibliography  on  the  subject,  with  more  than  350  references,  has  been 
collected  by  Dawson  and  Charlton  (1988).  Brief  review  of  some  of  the  previous 
work  follows. 

Pioneering  work  on  open-channel  flow  through  vegetation  was  performed  by 
Ree  and  Palmer  (1949;  see  also  Palmer  1945),  who  developed  a  method  for 
estimating  water  discharge  capacity.  They  employed  the  often-used  Manning’s 
coefficient,  concluding  that  the  n-URh  relationship  depends  on  the  physical 
properties  of  the  grass  and  is  thus  independent  of  channel  geometry  and  flow 
conditions.  Here  U  is  the  mean  velocity  and  Rh  is  the  hydraulic  radius.^ 

A  series  of  studies  has  been  conducted  at  the  University  of  Waterloo,  Canada, 
to  determine  the  flow  characteristics  of  vegetated  open  channels.  In  their  early 
investigation,  Kowen,  Unny,  and  Hill  (1969)  used  artificial  styrene  made 
roughness  elements  glued  to  the  bottom  of  a  laboratory  flume  to  study  flow  over 
simulated,  flexible  vegetation.  Pitot  tube  technique  allowed  them  to  measure 

*•  For  convenience,  symbols  and  unusual  abbreviations  are  listed  and  defined  in  the  Notation 
(Appendix  A). 
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velocity  distributions.  These  authors  found  a  good  fit  to  their  experimental 
results  with  a  modified  form  of  the  logarithmic  law  by  adjusting  the  values  of  the 
roughness  parameter  and  the  origin  intercept.  In  the  case  of  flexible  elements 
they  found  the  roughness  parameter  to  be  12  percent  larger  than  the  deflected 
height  of  the  elements.  They  suggested,  however,  the  use  of  a  log  law  using  the 
undeflected  height  of  the  plants,  and  by  comparing  it  with  experimental  field 
data  they  concluded  that  both  the  origin  intercept  and  the  slope  were  functions  of 
both  vegetation  density  and  plant  flexibility.  It  is  interesting  to  note  that  almost 
all  researchers  suggest  the  use  of  a  logarithmic  law  for  the  vertical  profile  of 
mean  velocities  above  the  plants,  hence  implicitly  assuming  the  existence  of  an 
equilibrium  layer,  i.e.,  with  production  of  turbulence  being  locally  balanced  by 
dissipation. 


Kowen  and  Unny  (1973)  conducted  a  series  of  experiments  simulating 
vegetation  by  using  plastic  strips  of  different  thicknesses.  They  proposed  the 
existence  of  three  basic  flow  regimes:  (a)  erect,  when  the  plastic  strips  are  erect 
and  stationary;  (b)  waving,  when  the  strips  undergo  a  waving  motion;  and  (c) 
prone,  when  the  strips  are  bent  over.  Similar  regimes  were  observed  by  Gourlay 
(1970)  for  Kikuyu  grass.  Despite  these  kinematical  classifications,  the  hydraulic 
behavior  of  simulated  vegetation  showed  the  existence  of  only  two  regimes, 
because  the  frictional  coefficient  for  both  erect  and  waving  “plants”  indicated 
identical  values,  whereas  much  lower  friction  factors  (by  a  factor  of  five)  were 
observed  for  the  prone  cases.  These  observations  clearly  indicate  the  existence 
of  a  common  turbulence  structure  for  flow  in  vegetated  channels  for  both  erect 
and  waving  plants,  whereas  a  different  turbulence  dynamic  probably  dictates  the 
behavior  of  prone  or  bent-over  vegetation.  This  latter  fact  is  herein  interpreted  as 
a  consequence  of  the  reduced  turbulent  diffusivity  coefficient  for  momentum  at 
the  top  of  the  plants  due  to  the  vertical  blockage  exerted  by  the  inclined 
elements.  These  investigators  also  introduced  a  stiffness  parameter,  MEI,  where 
M  is  the  relative  density  of  the  plants  and  El  is  the  stem  flexural  rigidity. 

Numerical  predictions  of  sediment  transport  capacities  in  vegetated 
free-surface  flows  were  attempted  by  Li  and  Shen  (1973)  based  on  a 
superposition  technique  for  the  wakes  generated  behind  isolated  elements,  a 
procedure  originally  proposed  by  Petryk  (1969).  They  assumed  local  drag 
coefficients  for  open  channels  of  about  1.2,  and  their  results  showed  mean  drag 
coefficients  close  to  1.1  for  staggered  arrangements  independent  of  plant  density, 
while  the  mean  drag  coefficient  for  square,  parallel  patterns  showed  increasing 
values  for  increasing  spacing.  They  applied  this  method  for  the  estimation  of 
bed  load,  and  compared  the  relative  effect  on  sediment  yields  by  various 
combinations  of  tall  vegetation. 
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Also  in  line  with  the  use  of  Manning’s  coefficient  as  a  measure  of  flow 
resistance,  Petryk  and  Bosmajian  (1975)  developed  a  quantitative  procedure  for 
predicting  this  coefficient  as  a  function  of  flow  depth  and  vegetation 
characteristics.  Their  method  considered  flow  depths  that  were  less  than  or  equal 
to  the  maximum  plant  height,  and  its  most  useful  application  is  in  predicting  the 
variation  of  Manning’s  n  with  depth. 

Among  the  attempts  to  build  models  that  provide  more  information  about  the 
flow  structure,  not  only  about  the  overall  flow  resistance,  Reid  and  Whitaker 
(1976)  developed  a  numerical  algorithm  for  wind-driven  flow  through  and  above 
vegetative  obstructions.  Accordingly,  they  divided  the  water  column  into  two 
layers,  one  within  the  canopy,  and  one  above  it  and  averaged  the  governing 
equations  within  each  layer.  The  main  drawback  of  this  approach  is  the  need  to 
specify  the  interfacial  stress  at  the  top  of  the  plants. 

In  the  area  of  atmospheric  boundary  layers,  Wilson  and  Shaw  (1977), 
recognizing  some  of  the  limitations  of  first-level  turbulence  closure  schemes, 
developed  a  higher-order  closure  model  for  atmospheric  flows  above  plant 
canopies.  At  the  same  time  these  authors  were  the  first  to  recognize  the 
necessity  of  spatial  as  well  as  temporal  averaging  of  the  governing  equations  for 
the  proper  one-dimensional  representation  of  the  problem. 

More  recently  Kowen  and  Li  (1980)  proposed  a  new  methodology  for  the 
design  of  channels  with  vegetative  linings,  thus  improving  the  traditional  n-URh 
method.  The  originality  of  this  new  method  consists  in  introducing  some 
biomechanical  concepts  and  proposing  a  field  methodology  for  estimating  the 
flexural  stiffness  of  natural  vegetation:  a  “board  drop  test”  and  a  vegetation 
height  method. 

Hino  (1981)  was  probably  one  of  the  first  researchers  to  address  the 
importance  of  vegetation  in  open  channels  from  an  ecohydrodynamic  point  of 
view.  He  also  pointed  out  the  particular  mathematical  and  numerical  difficulties 
that  arise  in  the  case  of  free-surface  flows,  which  transform  the  problem  into  a 
nonlinear  two-point  boundary- value  problem  with  an  implicitly  posed  boundary 
condition  at  the  bottom,  and  presented  some  numerical  results  as  well  as 
perturbation  solutions  using  a  mixing-length  closure  scheme. 

Raupach  and  Shaw  (1982),  based  on  previous  work  by  Wilson  and  Shaw 
(1977),  first  proposed  a  mathematical  procedure  for  obtaining  the  momentum 
and  energy  equations  in  multi-connected  flows  clearly  stating  the  rules  for  the 
commutation  of  spatial  averaging  operators  and  spatial  differentiation.  Their 
work  allows  for  the  identification  of  different  momentum  and  energy  dispersive 
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terms  arising  as  a  consequence  of  the  three-dimensional  nature  of  the  flow 
structure  as  well  as  of  the  noncommutation  of  the  operators  mentioned. 

One  of  the  few  laboratory  works  on  the  sedimentology  of  erect  vegetation  in 
open  channels  was  conducted  by  Tollner,  Barfield,  and  Hayes  (1982;  Tollner 
1974),  which  reported  good  predictions  of  sediment  transport  capacity  using 
parameters  similar  to  the  ones  proposed  by  Graf  (1971),  but  with  the  channel 
width  replaced  by  the  element  spacing.  Their  results  were,  however,  obtained  in 
a  relatively  short  (2.10  m)  and  narrow  (0.13  m)  channel,  where  the  achievement 
of  equilibrium  conditions  (at  least  for  suspended  sediment  profiles)  becomes 
questionable. 

The  primitive  nature  of  the  closure  scheme  used  by  Reid  and  Whitaker  (1976) 
motivated  Burke  and  Stolzenbach  (1983;  Burke  1982),  who  were  probably  one 
of  the  first  to  propose  the  use  of  a  two-equation  turbulence  closure  scheme  for 
free-surface  flows  through  obstructions.  In  this  kind  of  closure  the  eddy 
viscosity  is  assumed  proportional  to  the  product  of  a  characteristic  length  and  a 
velocity  scale,  both  obtained  by  solving  two  transport  equations,  so  that  they  do 
not  have  to  be  specified  a  priori  for  each  problem.  The  presence  of  vegetation 
was  accounted  for  in  the  turbulent  kinetic  energy  and  dissipation  equations  by 
the  introduction  of  drag-related  source  terms,  but  no  mathematical  derivation 
was  presented  to  justify  these  assumptions.  While  their  model  predictions  were 
generally  in  good  agreement  with  experimental  observations,  they  recognize  the 
lack  of  knowledge  about  the  value  of  the  drag  coefficient  of  the  elements  in 
open-channels  but  did  not  explain  satisfactorily  the  overestimation  of  the 
turbulent  kinetic  energy. 

Another  simpler  yet  useful  attempt  to  close  the  turbulence  problem  is  the  one 
due  to  Christensen  (1985),  who  used  the  mixing  length  approach  to  compute 
eddy  viscosities,  and  thus  developed  an  explicit  formula  for  the  velocity  profile 
over  a  flexible  roughness  layer  to  be  used  in  heavily  vegetated  rivers  and 
channels. 

In  light  of  the  proposed  averaging  procedure  by  Raupach  and  Shaw  (1982), 
Raupach  et  al.  (1986)  conducted  a  series  of  experiments  aimed  at  characterizing 
the  turbulence  structure  of  atmospheric  flows  over  a  vegetated  canopy.  They 
used  a  laboratory  wind  tunnel  with  a  model  plant  canopy  made  of  aluminium 
strips,  where  velocity  measurements  were  taken  at  several  points  both  above  and 
within  the  roughness  elements  using  a  special  three-dimensional  hot-wire 
anemometer.  They  were  thus  able  to  estimate  the  different  terms  composing  the 
turbulent  kinetic  energy  balance  within  the  canopy.  From  their  observations,  the 
importance  of  the  inertial  transport  term  atop  of  the  simulated  canopy  becomes 
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noticeable,  which  represents  a  major  loss  near  the  top  of  the  canopy  but 
constitutes  the  principal  gain  mechanism  lower  down. 

Saowapon  and  Kowen  (1989)  have  advanced  an  analytical  model  for 
predicting  vertical  velocity  profiles  in  vegetated  channels  that  accounts  for  the 
flexibility  of  the  plants.  While  the  results  of  both  the  Christensen  (1985)  and 
Saowapon  and  Kowen  (1989)  models  look  rather  encouraging  when  compared 
against  laboratory  observations,  it  is  clear  that  the  algebraic  scheme  used  to 
compute  eddy  viscosities  (i.e.  a  mixing  length  approach)  provides  only  limited 
information  on  the  effect  of  roughness  elements  on  the  diffusion  of  momentum 
(and  eventually  sediment). 

Kadlec  (1990)  obtained  a  power  law  resistance  function  for  overland  flow 
over  Spartina  grass  in  terms  of  depth  and  friction  slope.  The  exponent  on  the 
depth  appears  to  describe  both  the  vertical  vegetation  stem  density  and  the 
bottom  elevation  distribution  and  takes  a  value  close  to  three.  This  fact  seems  to 
indicate  a  strong  depth-dependent  behavior  in  wetlands,  where  depth-time 
variations  are  strongly  regulated  by  a  condition  of  small  and  slowly  varying 
depths. 

The  Kanazawa  University  group  (Tsujimoto  et  al.  199P  and  1991*’; 

Tsujimoto  1993;  Shimizu  and  Tsujimoto  1993)  has  reported  several 
open-channel  mrbulence  measurements  in  the  presence  of  vegetation.  A  series  of 
experimental  as  well  as  numerical  studies  has  been  conducted  concerning  both 
rigid  and  flexible  emergent  and  nonemergent  vegetation.  Their  results  for  rigid 
vegetation  show  that  an  almost  uniform  mean  velocity  distribution  prevails  when 
the  mean  flow  depth  is  smaller  than  the  vegetation  height,  with  negligible 
turbulent  momentum  exchange  and  small  turbulent  intensities.  On  the  other 
hand,  shear-dominated  flows  seem  to  prevail  for  nonemergent  vegetation,  even 
below  the  top  of  the  plants,  as  a  consequence  of  active  momentum  exchange 
between  the  faster  surface  flow  and  the  flow  within  the  simulated  vegetation.  In 
this  latter  case,  a  peak  in  the  Reynolds  stress  distribution  is  observed  at  the  top  of 
the  roughness  elements.  The  corresponding  results  for  flexible  vegetation  show 
that  the  mean  velocity  profile  is  no  longer  as  uniform  as  with  rigid  elements  for 
emergent  vegetation,  showing  slightly  decreasing  values  as  the  free  surface  is 
approached.  Turbulent  intensities  are  still  negligible.  On  the  other  hand, 
nonemergent  results  indicate  the  existence  of  a  deflection  point  near  the  top  of 
the  elements,  with  a  corresponding  peak  in  the  vertical  distribution  of  turbulent 
intensities.  Concerning  the  numerical  model,  relatively  good  agreement  with  the 
experimental  observations  was  obtained  with  a  modified  version  of  the  standard 
high-Reynolds-number  k-s  model,  also  with  drag-related  source  terms 
accounting  for  the  presence  of  vegetation.  Weighting  factors  in  these  source 
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terms  were  fit  to  reproduce  observed  distributions  of  mean  velocity  and 
Reynolds  stresses.  Although  the  numerical  code  is  very  similar  to  the  one 
employed  by  Burke  and  Stolzenbach  (1983),  the  weighting  factors  in  both 
schemes  are  radically  different,  the  most  striking  aspect  being  the  value  of  0.07 
obtained  by  Tsujimoto  et  al.  for  the  drag-related  turbulence  production  term 
compared  to  the  value  of  1.0  used  by  Burke  and  Stolzenbach.  It  is  worth 
mentioning  that  with  the  value  of  0.07,  Tsujimoto  et  al.  obtained  a  very  good  fit 
for  the  measured  turbulence  intensities.  In  a  more  recent  work,  however, 
Tsujimoto  and  Shimizu  (1994)  report  good  agreements  with  experimental 
observations  by  using  the  same  weighting  factors  as  Burke  and  Stolzenbach.  It 
is  therefore  not  clear  at  all  why  very  low  values  of  the  weighting  coefficients 
yield  good  agreement  with  the  observed  turbulent  kinetic  energy  profiles, 
whereas  from  a  physical  point  of  view  the  coefficient  of  the  work  done  by  the 
flow  against  form-drag  forces  should  be  equal  to  unity  (and  thus  very  close  to 
1.0  in  the  numerical  model).  Note  that  both  approaches  give  a  similar  degree  of 
fit  to  observed  profiles  of  mean  velocities  and  Reynolds  stresses,  which  should 
be  attributed  to  the  fact  that  eddy  viscosity  is  computed  as  proportional  to  the 
ratio  between  k  and  e,  so  that  certain  combinations  of  the  weighting  coefficients 
in  the  drag-related  source  terms  in  k  and  s  equations  yield  similar  values  of  eddy 
viscosity. 

Finally,  concerning  field  results.  Freeman,  Hall  and  Abraham  (1994) 
performed  several  field  tests  to  determine  Manning’s  n  values  and  sediment  trap 
efficiencies  for  stands  of  bulrashes.  Their  results  indicate  resistance  coefficients 
substantially  higher  than  the  ones  suggested  by  the  U.S.  Geological  Survey 
(Arcement  and  Schneider  1989).  Manning’s  n  was  observed  to  be  in  a  range 
between  0.26  and  0.70,  with  values  increasing  linearly  with  vegetation  density. 
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3  Theoretical  Considerations 


As  seen  in  the  previous  chapter,  most  attempts  to  numerically  simulate 
open-channel  flows  in  the  presence  of  vegetation  have  used  either  very 
simplified  or  more  complex  closure  schemes,  but  with  the  common  basis  of 
artificially  introducing  the  effect  of  vegetation  by  adding  body  forces.  As  will  be 
shown  in  this  chapter,  this  approach  has  led  to  some  inconsistencies  when 
numerical  results  are  compared  against  experimental  observations.  To  overcome 
such  difficulties,  the  governing  equations  are  herein  first  derived  for  a 
free-surface  flow  through  obstacles,  adapting  expressions  originally  developed 
for  atmospheric  flows  through  plant  canopies.  Through  this  derivation,  it  will  be 
shown  that  the  presence  of  vegetation  generates  dispersive  fluxes  of  momentum 
and  energy  as  well  as  viscous  and  form-drag  forces.  Since  the  latter  is  usually 
parameterized  using  a  drag  coefficient,  the  second  part  of  this  chapter  deals  with 
the  evaluation  of  this  coefficient  in  open  channels.  Once  the  governing 
equations  that  mathematically  define  the  problem  are  obtained,  the  third  part  of 
the  chapter  introduces  the  assumptions  needed  in  order  to  close  the  turbulence 
problem  at  first  level  using  a  two-equation  algorithm.  Finally,  some  limitations 
of  the  assumptions  introduced  are  presented,  followed  by  a  brief  description  of 
the  numerical  method  employed  to  solve  the  resulting  system  of  nonlinear  partial 
differential  equations. 

Governing  Equations  for  Fiow  Through  Vegetation 

This  section  deals  with  open-channel  flow  in  the  presence  of  vegetation. 
Therefore  the  analysis  presented  is  a  slight  modification  of  the  approach 
proposed  by  Raupach  and  Shaw  (1982)  for  atmospheric  flows  through  plant 
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canopies.  From  a  mathematical  point  of  view,  the  flow  of  water  through  and 
above  plants  presents  new  challenges  due  to  the  three-dimensionality  of  the 
turbulence,  thus  representing  a  highly  nonhomogeneous  flow  field.  Since  from 
an  engineering  perspective  a  one-dimensional  description  of  the  problem  is 
commonly  desirable,  the  need  for  spatially  averaging  (at  least  horizontally) 
naturally  arises  in  the  problem.  In  an  earlier  work,  Wilson  and  Shaw  (1977) 
already  noted  that  the  traditional  approach  of  arbitrarily  introducing  form-drag  as 
an  extra  body  force  in  the  momentum  equation  incorrectly  describes  the  effect  of 
wake  turbulence.  But  it  was  not  until  the  work  of  Raupach  and  Shaw  (1982)  that 
the  complete  set  of  equations  became  available. 

Wilson  and  Shaw  (1977)  offered  two  averaging  schemes  for  the  conservation 
equations.  In  the  first  one,  hereafter  termed  Scheme  I,  the  equations  describing 
the  instantaneous  flow  field  are  averaged  over  a  plane  large  enough  to  eliminate 
fluctuations  due  to  both  the  turbulent  scales  and  the  canopy  structure.  Consider 
an  open-channel  flow  through  a  regular  array  of  vertical,  rigid  cylinders 
simulating  vegetation  such  as  depicted  in  Figure  3.1.  Freezing  the  flow  at  any 
given  instant  and  analyzing  the  instantaneous  velocity  field  will  show  the 
existence  of  spatial  variations  due  to  the  canopy  (for  example,  differences 
between  u\,  Mj  and  Mj,  where  «( represents  the  component  along  the  i-axis,  x,-, 

of  the  instantaneous  velocity  at  location  j)  as  well  as  spatial  variations  at  similar 
locations  in  the  flow  due  to  the  intrinsic  nature  of  turbulence  (for  example, 
differences  between  u\,  u\  and  «]).  Therefore,  if  the  channel  is  sufficiently 
wide  and  long,  a  large  enough  horizontal  area  has  to  be  chosen  so  that  averages 
of  the  instjuitaneous  flow  field  performed  over  that  plane  will  provide  mean 
values  independent  of  spatial  variations  due  to  the  canopy  structure  and  the 
turbulence.  In  the  second  averaging  procedure  proposed  by  Wilson  and  Shaw 
(1977),  hereafter  called  Scheme  11,  the  three-dimensional  flow  structure  is  locally 
time-averaged  first  in  the  usual  way  to  filter  fluctuations  due  to  the  turbulence, 
and  then  spatially  averaged  to  eliminate  variations  due  to  the  canopy  structure. 
Referring  to  the  experiment,  at  each  location  there  would  be  a  fluctuating 
velocity  series  in  the  time  domain.  First  a  temporal  average  would  thus  be 
performed  at  each  location  to  get  rid  of  turbulent  fluctuations,  and  then  the 
time-mean  data  would  be  spatially  averaged  to  filter  spatial  variations.  It  is 
readily  seen  that  for  regularly  arranged  plants  the  extent  of  the  spatial  filter  in  the 
first  scheme  has  to  be  much  larger  than  in  the  second  one. 

In  the  following  discussion,  angle  brackets  and  overbars  will  indicate 
horizontal  and  temporal  averages,  respectively,  and  double  and  single  primes  will 
indicate  spatial  and  temporal  fluctuations,  respectively,  from  their  corresponding 
mean  values.  The  formal  definition  of  a  horizontal  average  of  a  variable  ^  is 
(see  Figure  3.1); 
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Figure  3.1  Schematic  of  flow  through  vegetation 


<  xp  > 
where 


xp(x,y)  dx  dy 


A  =  horizontal  area  of  rectangle  shown  in  Figure  3.1 


(1) 


From  a  mathematical  point  of  view  the  spatial  averaging  operator  in  a  flow 
through  obstacles  satisfies  all  but  one  of  the  commutation  properties  required  for 
a  turbulence  averaging  operator  (Monin  and  Yaglom  1971;  Schlichting  1979). 
This  exception  concerns  the  commutation  between  averaging  and  differentiation 
operations.  Raupach  and  Shaw  (1982)  clearly  show  how  Green’s  theorem  may 
be  used  to  demonstrate  that  if  W  is  constant  at  the  fluid-element  interface,  then 
horizontal  averaging  and  spatial  differentiation  commute  (i.e. 

<  dxp/dxi  >  =  d  <  xp  >  /dx,  ),  and  otherwise  they  do  not  commute.  In 
particular,  three  cases  of  interest  result:  spatial  differentiation  and  horizontal 
averaging  do  not  commute  for  pressure,  and  neither  do  Laplacian  operators  and 
horizontal  averaging  for  velocity,  but  first-order  spatial  differentiation  and 
horizontal  averaging  do  commute  for  velocity  and  its  higher  order  moments. 

Continuity  Equation 

The  instantaneous  continuity  equation  for  an  incompressible,  homogeneous 
and  steady  flow  is  (tensor  notation  will  be  used): 


^  =  0  (2) 

dXi 

Because  of  the  rules  enumerated  in  the  preceding  section,  both  schemes  yield 
essentially  identical  results  when  Equation  2  is  considered,  i.e.: 
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Scheme  I 


Scheme  11 


-r—  <  U:  >  =0 
dX:  ‘ 


4-  <^>  =0 
V.  * 


Momentum  Equations 


The  Navier-Stokes  equation  for  an  incompressible,  homogeneous  flow  is 


du:  du:  1  dp 

IF  q'^- 

where 
t  =  time 


p  =  instantaneous  pressure 


gi  =  component  in  the  ith  direction  of  the  gravitational  acceleration 
V  =  fluid  kinematic  viscosity 

The  resulting  momentum  equation  under  Scheme  I  is  found  by  following  the 
usual  way  of  first  replacing  >  +  u/  and  p  —  <  p  >  +  p",  and 

then  averaging  spatially  using  the  aforementioned  rules,  yielding  (Raupach  and 
Shaw  1982): 


d  <  U:  >  d  <  U:  >  Pi 

- ^ -  +  <  My  >  - r— ^ -  +  ^  <  My"  My"  >  = 

at  J  dxj  dXj  ^  J 


\  d  <  p  >  __  1  fdp^' 
Q  dxi  Q  \  dxi 


+  gy  +  V  <  My  >  +  V  <  V^My"  > 


Note  that  since  p"  is  not  constant  at  the  fluid-element  interface,  then 
<  dp''  jdx^  >  is  not  equal  to  d  <  p"  >  /8xy  (which  by  definition  is  zero). 

To  obtain  the  averaged  equation  under  Scheme  II,  first  time-average  the 
equation  in  the  usual  way,  then  substitute  mJ  =  <  m"  >  +  m^  and 
p  =  <  p  >  +  and  finally  a  spatial  average  is  performed  yielding: 


d  <  U:  >  _  d  <  U:  >  ^  ^  — - - 

— a7^  +  sjf-  = 


i  " e\^)  +  *.  +  ^  V  < «, >  +  V  <  > 
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It  is  easily  observed  that  Equations  5  and  6  are  identical,  provided  that  the 
averaging  plane  under  Scheme  I  is  large  enough  to  assure  that 

<  Ui  >  =  <  u~  >  and  <  p  >  =  <  p  > ,  with  the  only  exception  being  in  the 
form  of  the  Reynolds  stress  term.  Thus  the  one-dimensional  momentum 
equation  for  flow  through  obstacles  accounts  not  only  for  a  Reynolds  stress  due 
to  turbulent  momentum  transfer,  <  u-u-  > ,  but  also  for  stresses  that  arise  due 

to  spatial  variations  of  the  mean  flow  field,  <  ul  Wj  > .  Hence  the  total 
resulting  stress  becomes: 

<  Ui  Uj  >  =  <  Wi  TTj  >  4-  <  M,'  Uj  >  ^ ' 

Unfortunately  the  dispersive  fluxes  (first  term  on  the  right)  have  so  far  eluded 
direct  measurements  either  in  open  channels  or  atmospheric  boundary  layers,  so 
that  their  relative  effect  upon  the  total  stress,  albeit  believed  to  be  small,  still 
remains  unknown.  From  a  mathematical/physical  point  of  view  it  becomes 
however  clear  that  the  simple  addition  of  drag-related  body  forces  in  the 
momentum  equation  is  essentially  incorrect  since  the  dispersive  fluxes  are  not 
included.  The  problem  becomes  even  more  relevant  when  higher-order  moments 
of  velocity  are  considered. 

Energy  (Second-Order  Moment)  Equations 

The  usual  procedure  for  obtaining  the  equations  for  the  mean  flow  as  well  as 
for  the  turbulent  kinetic  energy  (e.g.,  Hinze  1975)  will  be  followed.  Thus,  the 
total  kinetic  energy  under  Scheme  I  yields: 


^  > 


U:  >  <  U:  > 


+  2  ^ 


U;  U; 


(8) 


With  due  regard  to  the  commutative  properties  mentioned  in  the  previous 
section,  the  equation  for  the  mean  kinetic  energy  is  (Raupach  and  Shaw  1982): 


a  a  \  <  M;  >  <  M,-  > 

-h  <  M;  >  ®  *  '  ' - 


dt 


d  <  U:  > 

<  >  - r-^ - 

^  J  dXi 


dXj^ 


<  Ul  >  <  u/'  Uj"  >  + 


p  >  <  Uj  > 


+  v  <Ui><  i  <  «,■  > 


(9) 
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and  likewise  the  budget  of  turbulent  kinetic  energy  under  such  scheme  gives 
(Raupach  and  Shaw  1982): 
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3?  +  ^  ex, 

<  J 


d  \<  Ui  Ui  > 


=  -  <  U^Uj  > 


d  <  Ui  > 
dxj 


3  /<  >  <P"^J 


> 


dX: 


+ 


+  V  <  > 


+  77  <  M;  > 


(10) 


The  right-hand  side  of  each  equation  contains  four  terms:  (a)  a  shear-production 
term,  a  source  term  in  Equation  10  but  a  sink  term  in  the  budget  of  mean  flow 
energy,  which  converts  mean  kinetic  energy  to  large-scale  turbulent  kinetic 
energy;  (b)  a  turbulent  transport  term  with  the  usual  inertial  and  pressure 
components;  (c)  a  viscous  term;  and  (d)  a  wake-production  term,  a  source  term  in 
Equation  10  but  a  sink  in  Equation  9,  representing  the  rate  of  work  of  the  mean 
flow  against  the  force  exerted  by  the  obstacles.  The  viscous  term  accounts  for 
molecular  diffusion,  molecular  transport,  and  viscous  dissipation  of  turbulent 
kinetic  energy  (Townsend  1976).  The  fourth  term  in  Equation  10  accounts  for 
the  conversion  of  both  mean  and  large-scale  turbulent  kinetic  energy  toward 
smaller  scale  turbulent  kinetic  energy  in  the  wakes  of  the  elements,  which  is 
sometimes  referred  to  as  “short-circuited  cascade”  (Raupach  and  Thom  1981). 
This  wake-generated  turbulent  kinetic  energy  has  therefore  a  scale  proportional 
to  the  dimensions  of  the  elements  in  the  canopy,  i.e.,  much  smaller  than  the 
typical  length  scales  of  shear-generated  eddies  (Raupach  and  Thom  1981; 
Raupach  and  Shaw  1982). 

Under  averaging  Scheme  n,  the  decomposition  of  the  total  kinetic  energy  is  a 
little  more  complicated,  yielding: 


i  i  <  M,.  «,•  >  +  ^  <  M,'m/  > 


(11) 


1 


1 


uj>  <  u:>  +  IT^  iq  >  +  ^<  > 


where  the  last  two  terms  on  the  right-hand  side  represent  different  components  of 
the  turbulent  kinetic  energy  under  this  scheme. 


Budgets  for  each  of  the  last  two  terms  on  the  right  in  Equation  1 1  are  readily 
obtained  (Raupach  and  Shaw  1982): 
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—  <  uluj  > 


^  d  a  \  > 

—  +  <  zr  >  - - — ’ — 

,dt  J  dX:  2 


d  <  Ui> 
dXj 


^  /  <  u/Ui  U-  >  <  u/u/  Uj  >  <  p'uj  > 

2  2  ■*■  Q 


+  V  <  m/VV  >  - 

^  ^  ^  J  dXj  I 

and 


a  a  \  ^  ^ 

“  4*  <  zr  >  —  — - — ’ — 
dt^  ^J  dX:  I  2 


d  <  U:  > 
—  <  W;  M;  >  - T - 

'  J  dX: 


+  (  M;  Uj  -T - 


dXj 


<  U-  U:  U:  >  <  U:  U:  U:  >  <  f  U:  > 

- r - 1 - T - 1 - T; - 


+  V  <  U:  V^w/'  >  +  77  <  M,-  >  y  » 

^  *  Q  '  \  a;c; 


dp' 


(12) 


(13) 


The  four  terms  on  the  right-hand  side  of  Equation  12  have  similar  meanings 
to  the  ones  in  Equation  10,  except  that  the  wake-production  term  appears  here  as 
a  horizontal  average  of  the  product  of  local  deviations  of  Reynolds  stresses  and 
velocity  gradients  from  their  spatial-averaged  values.  On  a  smaller  scale,  this 
last  term  produces  turbulent  kinetic  energy  in  the  same  way  as  does  the 
shear-production  term.  More  attention  has  to  be  given  to  the  five  terms  on  the 
right  in  Equation  13:  (a)  a  production  term;  (b)  the  wake-production  term  of 
Equation  12,  here  a  sink  term;  (c)  a  turbulent  transport  term  involving  the  role  of 
dispersive  fluxes  of  energy;  (d)  a  viscous  term;  and  (e)  a  wake-production 
(source)  term,  similar  to  the  fourth  term  in  Equation  10. 

Careful  analysis  of  the  preceding  expressions  allows  for  better  insight  into  the 
turbulence  stmcture  and  its  generation  mechanisms  in  flows  through  vegetation. 
Basically  it  can  be  observed  that,  irrespective  of  the  averaging  scheme,  the 
budget  of  turbulent  kinetic  energy  is  composed  of  sources,  sinks,  and  transport 
terms.  Two  characteristic  processes  act  as  turbulent  kinetic  energy  generators, 
i.e.,  transferring  energy  from  larger  scales  (either  mean  flow  or  larger  eddies) 
toward  turbulent  fluctuations  in  space  or  time  at  smaller  scales:  (a)  the  work  of 
Reynolds  and  dispersive  stresses  against  mean  velocity  gradients;  and  (b)  the 
work  of  mean  flow  or  large  eddies  against  pressure  differences  due  to  the 
obstacles.  Looking  at  Equation  7  and  at  the  first  term  on  the  right  of  Equation 
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10,  we  can  see  that  the  action  of  mechanism  (a)  may  in  turn  be  subdivided  as 
d  <  TT;  > 


<  w/  > 


dxj 


which  is  the  same  as  the  first  term  on  the  right  of  Equation  12,  and  therefore 
contributes  to  the  generation  of  fluctuations  in  time,  and 
d  <  TT;  > 


<  U:  U:  > 


dxj 


which  is  equal  to  the  first  term  on  the  right  of  Equation  13,  and  thus  generates 
spatial  perturbations  of  time-averaged  values.  On  the  other  hand,  the  work  of  the 
mean  flow  against  pressure  differences  in  space  (i.e.,  mechanism  b)  is  a  source 
term  for  the  budget  of  spatial  fluctuations  of  time-averaged  velocities,  where  a 
shear-generation-like  term  appears  as  a  sink,  thus  transferring  energy  from  space 
fluctuations  toward  small-scale  fluctuations  in  time. 


Regarding  transport  processes,  the  second  term  on  the  right  in  Equation  12  is 
identical  to  the  corresponding  term  in  the  turbulent  kinetic  energy  budget  of 
shear  flows  without  obstacles,  with  the  exception  being  the  appearance  of  a 

dispersive  flux  of  turbulent  kinetic  energy,  <  u/  > . 


In  the  turbulent  kinetic  energy  budgets  there  are  two  viscous-related  sink 
terms,  accounting  for  the  direct  conversion  of  mechanical  energy  into  heat.  The 
one  in  Equation  12  is  related  to  the  spatial  average  of  the  typical  viscous 


du/  ( du/ 

dissipation  of  turbulent  kinetic  energy,  e  =  I  1 

This  relation  can  be  shown  as  follows  (Townsend  1976;  Hinze  1975): 


,  d^u/ 


9jc? 

j 


\ '2] 

d^U:'  u/ 

-J-  _  J 

1 - 

'  dx^dxj 

—  a 


(14) 


so  that  at  high  enough  Reynolds  numbers  the  viscous  term  in  Equation  12  is 
equivalent  to  the  rate  of  dissipation  of  turbulent  kinetic  energy  into  heat,  hence 
determining  the  viscous  cutoff  of  turbulent  fluctuations  in  time.  The  other 
viscous  term,  the  fourth  term  in  Equation  13,  accounts  for  the  direct  dissipation 
into  heat  of  spatial  fluctuations  of  time-averaged  mean  velocities. 

There  are  two  limiting  cases  worth  being  analyzed.  The  first  one  is 
considered  in  the  work  of  Raupach  and  Shaw  (1982)  and  concerns  the  case  when 
the  length  scale  of  the  canopy  elements  (and  of  their  wakes,  or  in  other  words  the 
scale  of  the  wake-generated  turbulence)  is  much  larger  than  the  Kolmogorov 
microscale,  ?/,  so  that  the  viscous  term  in  Equation  13  becomes  negligible.  In 
this  situation,  if  all  the  dispersive  fluxes  are  considered  to  be  of  lower  order  of 
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magnitude,  then  for  steady  advection-free  conditions: 


> 


(15) 


In  other  words,  the  work  of  the  mean  flow  against  pressure  differences  becomes 
equal  to  the  wake-production  term  for  the  turbulent  fluctuations  in  time. 


The  other  limiting  case  is  when  the  length  scale  of  the  canopy  elements  (and 
of  their  wakes,  i.e.  the  scale  of  the  wake-generated  turbulence)  is  much  smaller 
than  (or  even  of  the  order  of)  the  Kolmogorov  microscale.  In  this  situation 
almost  all  the  energy  arising  from  the  work  of  the  mean  flow  against  pressure 
forces  is  spent  in  the  generation  of  spatial  fluctuations,  and  is  therefore  directly 
dissipated  into  heat.  In  steady  advection-ffee  conditions,  it  follows  that: 


—  V  <  M,  > 


> 


(16) 


So  that 


Uj  Uj 


dXj 


0 


and  hence  there  is  a  negligible  contribution  from  the  wakes  to  the  spatial  average 
of  the  turbulent  fluctuations  in  time.  The  first  of  these  two  situations  seems  to 
be  common  to  atmospheric  flows,  whereas  the  second  situation  is  more  common 
to  water  flows  with  relatively  low  plant  concentrations.  This  is  reasonable, 
considering  that  the  Kolmogorov  microscale  is  smaller  in  air  than  in  water.  In 
addition,  the  characteristic  length  scales  of  canopy  elements  in  atmospheric 
flows  can  be  expected  to  be  in  general  much  larger  than  those  found  in  water 
flows. 


The  discussion  in  the  previous  paragraphs  clarifies  the  problem  mentioned  in 
the  literature  review  concerning  the  different  coefficients  assigned  to  the 
wake-production  terms  in  different  turbulence  models.  It  becomes  therefore 
obvious  that  if  one  is  trying  to  numerically  simulate  the  spatial  average  of  the 

local,  time-averaged  turbulent  kinetic  energy  (or  any  <  >  for  that  matter) 

in  a  flow  with  elements  of  the  order  of  the  Kolmogorov  microscale,  then  the 
wake-related  source  term  in  the  energy  equation  would  be  almost  zero.  In  other 
words,  in  this  case  the  drag-related  weighting  factors  in  the  turbulent  kinetic 
energy  and  in  the  dissipation  equations  would  be  very  close  to  zero.  However, 
for  the  numerical  computation  of  the  total  turbulent  kinetic  energy  (i.e. 

(<  IT-  uj"  >  -f  <  «/  «/  >)/2),  these  coefficients  are  expected  to  be  close  to 
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1.0  and  1.33,  respectively  (see  “Turbulence  Modeling  by  a  First-Level 
Two-Equation  Closure  Scheme”). 

In  the  computation  of  flow  through  vegetation  using  a  two-equation  model,  a 
question  arises  about  which  turbulent  kinetic  energy  will  be  simulated:  the  one 
governed  by  Equation  10  or  12?  Moreover,  since  in  two-equation  models  the 
transport  equation  for  turbulent  kinetic  energy  provides  a  velocity  scale  for 
estimating  the  turbulent  eddy  viscosity,  which  one  of  the  two  will  provide  a 
better  approach  for  this  purpose?  The  answer  to  this  question  will  also  determine 
the  type  of  transport  equation  used  for  estimating  the  viscous  dissipation  rate  of 
turbulent  kinetic  energy.  In  other  words,  which  will  be  simulated 

V  <  >  or  V  <  >  ?  But  before  these  questions  are  answered, 

i.e.,  before  addressing  the  point  on  the  modeling  of  the  former  expressions,  some 
considerations  concerning  the  parameterization  of  drag  forces  will  be  presented. 

Modeling  of  Form-Drag  Forces 
in  Open-Channel  Flows 

The  previous  section  discussed  how  drag-related  terms  can  be  introduced  in 
the  conservation  equations  without  arbitrarily  introducing  body  forces.  The 
present  section  will  deal  with  the  modeling  of  such  forces.  As  mentioned  before, 

the  term  ^{dp" / dx^  represents  the  so-called  drag  force  per  unit  volume.  To 

demonstrate  this  assertion  Figure  3.2  shows  a  schematic  of  the  pressure  field  for 
an  isolated  two-dimensional  object.  By  definition  the  pressure  force  (per  unit 
length  in  z),  fx,  on  the  perimeter,  s,  of  the  cylinder  acting  in  the  ^-direction  is 
(e.g.  Pantom  1984): 

fx  =  i  nxp  ds  (17) 


where 

n  =  vector  normal  to  the  perimeter  of  the  object 
%  =  ;c-component  of  vector  n 

It  is  readily  seen  that  at  a  fixed  spanwise  location,  the  longitudinal  gradient  in  p" 
times  Axis  equal  to  (tixuP  +  Hxd  p),  where  and  tixd  represent  the  vector  Hx  in 
the  upstream  and  downstream  faces  of  the  object,  respectively,  and^x  the 
distance  between  these  two  points  in  the  x-direction. 

In  fluid  mechanics,  the  drag  force  is  usually  parameterized  as; 


i 

Q 


-  ^  a 


<  M,-  >  2 


(18) 
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Figure  3.2  Schematic  of  pressure  field 


where  a  is  the  ratio  between  the  sum  of  the  differential  frontal  areas  of  the 
obstacles  divided  by  the  differential  volume  of  fluid  (Figure  3.3)  and  thus  has 
dimensions  of  L"-^,  and  Qj  is  the  so-called  drag  coefficient  which  physically  is 
proportional  to  the  momentum  thickness  of  the  wake  behind  the  object. 

Although  the  determination  of  Qj  is  a  key  factor  in  the  modeling  of  flow 
through  obstacles,  very  few  experimental  observations  exist  concerning  the 
determination  of  this  coefficient  in  the  particular  case  of  open-channel  flows. 
Realizing  this  problem,  during  the  completion  of  the  present  work  a  set  of 
experiments  was  conducted  at  the  Hydrosystems  Laboratory,  University  of 
Illinois  at  Urbana-Champaign,  in  order  to  specify  values  of  Qj  for  free-surface 
flows  through  simulated  vegetation  (Dunn  1996).  Rigid  as  well  as  flexible 
cylinders  were  used  in  the  study.  Using  a  new  methodology  developed  to 
evaluate  the  drag  coefficient  based  on  vertical  profiles  of  spatial  and  temporal 
mean  velocity  and  Reynolds  stresses,  results  showed  that  Cu  =  1.12)  ±  15 
percent,  for  the  range  of  dimensionless  parameters  employed  in  the  work. 


4 

Figure  3.3  Definition  diagram  for  cylinder  density 
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Experimental  data  from  this  study  will  be  compared  herein  with  results  from  the 
numerical  simulations. 


Turbulence  Modeling  by  a  First-Level, 
Two-Equation  Closure  Scheme 

As  can  be  observed  from  previous  expressions,  transport  equations  for 
velocity  moments  of  any  order  involve  knowledge  of  higher  order  moments,  in  a 
way  that  produces  a  mathematical  problem  with  more  unknowns  than  equations. 
This  is  the  so-called  closure  problem  in  turbulence. 

In  general  three  different  alternatives  are  available  today  for  the  numerical 
computation  of  turbulent  flows: 

a.  Direct  numerical  simulation  using  the  full  set  of  Navier-Stokes  equations,  hence 

simulating  all  eddy  sizes  (e.g.,  Kim,  Moin  and  Moser  1987). 

b.  Explicit  simulation  of  only  the  large,  energy-containing  eddies,  which  have 
length  scales  determined  by  each  particular  problem  (resolved  scales),  with  the 
flux  of  energy  towards  the  smaller  eddies  in  the  spectrum  (subgrid  scales) 
modeled  by  introducing  an  effective  viscosity,  which  increases  the  molecular 
viscosity  of  the  fluid  (e.g.  Piomelli  1994). 

c.  Use  of  the  Reynolds-averaged  form  of  the  Navier-Stokes  equations  (or  similar 

equations  for  higher  order  moments)  plus  some  assumptions  that  allow  solving 
the  closure  problem  of  having  more  unknowns  than  equations  (Rodi  1984). 


The  closure  problem  at  the  first-order  level,  that  is  in  the  turbulence-averaged 
form  of  the  Navier-Stokes  equations,  has  been  traditionally  solved  by  means  of 
different  numbers  of  equations  having  as  a  common  framework  an  eddy  viscosity 
model.  Thus  this  approach  inevitably  breaks  down  where  the  concepts 
underlying  the  eddy  viscosity  hypothesis  are  in  violation  of  the  physical 
processes  (see  section  in  this  Chapter  ’’Limitations  of  Turbulence  Models  Based 
on  Flux-Gradient  Approximations).  Since  the  late  seventies,  models  have  also 
been  developed  in  engineering  practice  for  higher  order  closures.  Second-order 
models,  for  example,  are  based  on  the  full  equations  for  the  Reynolds  stress 
tensor,  and  of  course  third-order  moments  are  being  modeled  (Wilson  and  Shaw 
1977). 

According  to  the  number  of  transport  equations  being  used  for  closing  the 
problem  of  the  Reynolds  stresses,  the  models  have  been  called  zero-,  one-  and 
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two-equation  models.  Since  all  these  models  rely  on  the  concept  of  an  eddy 
viscosity,  and  its  determination  requires  the  knowledge  of  one  velocity  and  one 
length  scale,  the  problem  reduces  to  the  estimation  of  these  two  variables  in  the 
turbulence  field.  The  simplest  prescription  of  the  Reynolds  stress  in  the  level  of 
zero  equation  (or  algebraic  models)  is  the  well  known  mixing-length  model, 
obtained  by  applying  the  methods  of  gas  kinetic  theory  to  turbulent,  macroscopic 
motions  of  fluid  continuum  (e.g.,  McComb  1990).  Although  successfully 
applied  in  many  situations,  its  major  drawback  is  in  its  lack  of  universality:  i.e., 
the  prescription  of  the  mixing-length  varies  from  one  type  of  flow  to  another. 

The  next  order  of  difficulty  is  the  one-equation  model,  which  makes  use  of  a 
transport  equation  for  the  turbulent  kinetic  energy  (assumed  proportional  to  the 
square  of  the  characteristic  velocity)  together  with  some  assumptions  to  prescribe 
the  production,  diffusion,  and  dissipation  terms.  In  this  case  a  length  scale  still 
needs  to  be  specified  by  means  of  some  empirical  relation. 

The  next  level  of  complication  is  the  introduction  of  a  second  transport 
equation  with  the  help  of  which  the  required  length  scale  can  be  computed.  The 
basis  for  these  kinds  of  models  seems  to  have  been  given  some  fifty  years  ago 
almost  simultaneously  by  Kolmogorov  (1942)  and  Prandtl  (1945).  According  to 
Kolmogorov  (Barenblatt  1995)  at  any  point  of  a  turbulent  flow  the  statistical 
dimensionless  properties  of  the  vortex  dissipative  stmctures  are  similar,  and  only 
their  time  and  length  scales  are  different.  Both  scales  may  be  estimated  by 
different  sets  of  two  transport  equations,  either  for  (k,s),  (Jc,o}),  etc.  (k,s,l,co 

representing  the  turbulent  kinetic  energy,  dissipation  rate,  length  scale  and 
dissipation  per  unit  turbulent  kinetic  energy,  respectively).  The  first  of  the 
former  set,  the  k-e  model,  is  probably  the  most  commonly  used  model  in 
engineering  practice,  and  has  proved  to  be  a  reliable  tool  in  a  wide  variety  of 
problems  in  hydraulic  and  environmental  engineering  (Rodi  1984). 

In  the  k-£  model  the  Reynolds  stresses  are  estimated  using  the  eddy  viscosity 
concept  as  follows: 


V  Y 


ydXj  dXij 


3  ^ 


where 

vj  =  kinematic  eddy  viscosity 


V  Y 


Cfi  =  parameter  with  standard  value  of  0.09 


5;-  =  Kronecker  delta 

A  similar  approach  will  be  followed  herein,  namely: 


(19) 
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(20) 


—  <  Ui  Uj  >  =  —  <  M;  Uj  >  —  <  «/  Uj  > 


I  d  <  U:  >  d  <  Uj  > 


-  h  * 


The  question  then  concerns  the  turbulent  kinetic  energy  to  be  modeled  (i.e., 
k  =  <  Uj  Uj  >  jl  or  k  =  <  u'i'u'i  >  /2 ). 

The  answer  to  this  question  can  be  obtained  by  contracting  Equation  20,  i.e., 
making  i=j,  yielding: 


<  Uj 


> 


<  U:  U:  >  = 


(21) 


From  continuity  (Equation  3)  the  first  term  on  the  right  vanishes,  so  that 
Equation  21  reduces  to: 


,  _  <  Uj  Uj"  >  +  <  Uj'  Uj'  >  (22) 

^  “  2 

This  finding  clarifies  which  rate  of  turbulent  dissipation  has  to  be  modeled. 
Since  by  definition  the  characteristic  velocity  scale  is  considered  proportional  to 
the  turbulent  kinetic  energy  defined  in  Equation  22,  then  the  dissipation  rate 
defined  in  Equation  10  has  to  be  used  accordingly  for  defining  the  associated 
length  scale,  namely  £  =  v  <  M;  V^m/’  >. 

One  last  consideration  is  needed  in  order  to  model  Equation  10.  The  exact 
form  of  the  inertial  and  pressure  transport  terms  is  of  no  practical  use  since  it 
involves  unknown  correlations  of  higher  order.  To  obtain  a  closed  set  of 
equations,  assumptions  similar  to  those  used  in  the  standard  k-E  model  are  made, 
namely  that  the  total  (inertial  and  pressure)  diffusive  flux  of  k  can  be  assumed 
proportional  to  its  gradient: 


<  Uj"uj"uj"  >  <  p"uj"  >  Vj-  d  <  Uj"Uj"  >  /2  ^23 

2  2  ""  ^  ^ 

Now  the  set  of  partial  differential  equations  that  will  be  numerically  modeled 
for  simulating  the  turbulence  structure  of  uniform,  two-dimensional 
open-channel  flow  through  vegetation  can  be  written  down  as  follows: 


a.  Continuity  Equation: 


d  <  Uj  > 

dx 


=  0 


(24) 
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b.  x-Momentum  Equation: 
d  < 

dt 

where 

g  =  gravitational  acceleration 


u,  >  a  d  <  u,  > 

Tt -  =  ^^0  +  +  — Yz — 


(25) 


So  =  bed  slope 

c.  Turbulent  Kinetic  Energy  Equation: 


dt  dz  [Vta  ^  dz_ 


+  -  S  +  Cf^fx  <  Ui  > 


(26) 


d.  Dissipation  Rate  Equation: 


de 

dt 


dz 


(^  +  V)^ 


(27) 


where: 


and  the  set  of  standard  constants  takes  the  following  values:  =  0.09,  Cj  = 
1.44,  C2  =  1.92,  oit  =  1.0  and  =  1.30.  The  parameters  and  have  to  be 
modeled  and  are  sometimes  considered  results  of  the  model  calibration 
(Tsujimoto,  Kitamura  and  Okada,  1991^).  The  unsteady  terms  in  the  previous 
equations  are  retained  only  for  computational  purposes,  so  that  a  steady  solution 
is  reached  as  an  asymptotic  state  (see  the  section  in  this  Chapter,  “Numerical 
Algorithm”). 

However,  comparing  Equations  10  and  26,  one  expects  the  value  of  the 
coefficient  to  be  equal  (or  very  close)  to  one.  Moreover,  it  can  be  shown 
(Burke  1982)  that  for  the  s-equation  to  be  in  balance,  the  value  of  the  coefficient 
Cfi:  has  to  be  dependent  upon  the  value  of  C^.  To  clarify  this,  consider  a  steady, 
horizontal  flow  through  vertical,  infinite  long  cylinders,  where  all  derivatives  in 
the  vertical  direction  vanish.  Then,  from  the  k-equation  e  =  C^fx  <Uy  >, 

and  from  the  8-equation  Cj  C^fx  <Ui>  =  C2£sothat,  =  C^jC-^  C^. 
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In  order  to  solve  the  system  of  partial  differential  equations,  appropriate 
boundary  conditions  for  each  variable  have  to  be  specified.  This  issue  is 
addressed  in  the  next  section. 


Boundary  Conditions  and  Constants 
for  Open-Channei  Fiow  Through  Vegetation 


One  of  the  biggest  limitations  of  the  set  of  partial  differential  equations  given 
in  the  previous  section  is  that  viscous  effects  have  been  neglected,  and  thus  the 
model  is  not  able  to  resolve  flow  regions  too  close  to  solid  boundaries.  In  other 
words,  the  model  is  expected  to  yield  acceptable  results  only  in  local, 
high-Reynolds-numbers  conditions.  In  the  standard  version  of  a 
high-Reynolds-number  k-8  model,  values  of  velocity,  turbulent  kinetic  energy, 
and  dissipation  are  specified  at  a  point  near  the  wall,  located  in  the  so-called 
equilibrium  region,  where  the  flow  exhibits  such  large  local  rates  of  energy 
production  and  dissipation  that  both  terms  are  approximately  in  local 
equilibrium.  Basically  this  assumption  yields  values  for  velocity,  k  and  s  related 
to  the  existence  of  a  semi-logarithmic  mean  velocity  profile.  On  the  other  hand, 
the  free  surface  region  is  sometimes  treated  as  a  symmetry  plane  (i.e.,  the  fluxes 
of  all  variables  are  zero,  which  is  known  as  “rigid  lid  assumption”),  but  more 
accurate  results  are  obtained  when  turbulence  damping  effects  are  considered  by 
specifying  values  of  the  dissipation  rate  as  a  function  of  k  and  the  flow  depth  H 
(Celik  and  Rodi  1984, 1988).  This  latter  approach  is  in  line  with  experimental 
observations  that  show  e  to  be  proportional  to  the  ratio  Urms^^Lx>  where  Urms  is 
the  root-mean-square  value  of  the  streamwise  velocity  fluctuations  and  is  their 
macro-length  scale  (approximately  constant  and  equal  to  70  percent  of  the  flow 
depth  in  the  free-surface  region,  Nezu  and  Nakagawa  1993).  This  latter  approach 
will  be  followed  in  the  present  work,  so  that  the  boundary  conditions  to  be  used 
are: 


a.  At  the  bed: 


b.  At  the  free  surface: 


M  =  dk  ^  r.  _  (29) 

dzdz  b,H 


where 
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=  model  coefficient 


Zo  =  first  grid  point  away  from  the  wall 

E  =  roughness  parameter,  approximately  equal  to  9  for  hydraulically 
smooth  conditions  and  to  Wn>/m/ks  for  fully-rough  beds. 

ks  =  equivalent  sandgrain  roughness 

It  is  worth  mentioning  that  Equation  29  is  only  valid  under  non-emergent 
conditions.  Garcia  (1992)  showed  how  these  boundary  conditions  can  be 
modified  to  account  for  buoyancy  effects  induced  by  sediment  in  suspension. 

Concerning  the  value  of  the  constant  C^,  different  approaches  exist  in  the 
literature.  Celik  and  Rodi  (1984)  reported  that  values  of  (!^  =  0.05  result  in 
predictions  of  near-bed  values  of  streamwise  velocity  and  kinetic  energy  in  good 
agreement  with  experimental  observations.  Rodi  (1976)  proposed  an  algebraic 
expression  for  estimating  Reynolds  stresses,  from  which  a  formula  for  Cji  can  be 
obtained,  which  shows  this  value  to  be  a  function  of  the  ratio  between  production 
and  dissipation  of  turbulent  kinetic  energy.  In  modifying  this  expression  to 
account  for  wake-generated  turbulence,  a  new  expression  for  will  be  obtained 
in  the  next  section. 

Estimation  of  Reynolds  Stress  Tensor  Components 

As  mentioned  in  the  preceding  section,  Rodi  (1976)  proposed  an  algebraic 
expression  for  estimating  the  different  components  of  the  Reynolds  stress  tensor. 
His  approach  is  slightly  modified  in  what  follows  to  account  for  wake-generated 
turbulence.  Consider  the  transport  equation  of  k\ 


Dk 

Dt 

Dif{k)  +  P  —  e  -k  Pw 

where 

Dif(k) 

=  the  diffusive  transport  of  k 

D(.)/Dt 

=  total  derivative 

Pw 

=  wake  production  term 

A  similar  expression  can  be  obtained  for  each  Reynolds  stress: 


(30) 


D  <  U,"Uj"  > 

Dt 


Dif(<  >)  -I-  Py 


(31) 


where  the  third  and  fourth  terms  on  the  right  account  for  pressure-strain  effects 
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(Rodi  1976;  see  also  Launder,  Reece  and  Rodi  1975).  As  it  can  be  observed. 
Equations  30  and  3 1  are  differential  equations  due  to  their  left-hand  side  and  the 
first  term  on  the  right  of  each  expression.  Now,  mathematically  it  can  be  written: 


D  <  ui'uj  > 
Di 


Dt 


+  k- 


Dt 


(32) 


where  a  similar  expression  can  be  obtained  for  Dif( <Ui  uj  > ).  If  the  ratio 

ff  n 

<Ui  >A  can  be  assumed  approximately  constant  in  the  computational  domain 
(something  that  is  in  fairly  good  agreement  with  the  authors’  own  experimental 
observations,  as  it  will  be  shown  later),  it  is  therefore  possible  to  write: 


D  <  u/'uj"  > 
Dr 


-  Dif{<  uruj"  >) 


UiUj 

'Dk 

k 

[Dt 

Difik) 


and  hence  from  Equation  30  it  follows  that: 


(33) 


D  <  >  <  m/'m,"  > 

- ^ - Difi<  uru/'  >)  « - '-jf - (P  -e  +  Pw) 


(34) 


Combining  Equations  31  and  34  yields: 


<  u/'uj"  > 


iP-e  +  Pw)  =  P,J  -  C,4(<  uru/'  >  - 


(35) 


or: 


<  > 
k 

where 

njf  ^  i  j _ L  ( P  +  __  1 

^  ^  I  f  ^  /  (3' 

Equation  36  is  an  algebraic  expression  for  obtaining  the  Reynolds  stresses,  once 
P,  £,  Pij,  Pw  and  Pwij  are  known. 

An  important  consequence  of  Equation  36  is  that  if  the  value  of  <ui"u3''>  is 
computed,  and  assumed  to  be  equal  to: 


f'^y  ^  +  c,^(«  0 


1  -  7 
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€  3^ijs 
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PWr, 


(36) 
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-  <  m/'ms"  > 


(38) 


= 


5  <  Wj  > 

dz 


_  ^  yfcza  <  Mj  > 

-  dz 


then  it  is  easy  to  show  that: 


Cfi 


2  i-y  1  - ^(1  - yP/e) 

3  Cj^ 


(39) 


which  clearly  shows  to  vary  as  a  function  of  both  ratios  P/e  and  Pw/e. 


Regarding  the  values  of  the  coefficients  y  and  Qr  (see  Launder,  Reece  and 
Rodi  1975),  the  former  takes  a  value  of  0.60  for  isotropic  turbulence,  whereas 
the  latter  was  originally  found  to  be  equal  to  1.4  by  Rotta  (1951).  However, 
Rotta  (1962)  later  showed  that  a  value  about  twice  as  large  provided  a  better  fit 
to  Uberoi’s  (1957)  data  on  the  decay  of  highly  anisotropic  turbulence.  Rodi 
(1976)  suggests  the  use  of  y=0.4  and  Cir=2.5.  Figure  3.4  illustrates  the 
variation  of  Qt  with  the  ratio  between  total  production  and  dissipation,  for  the 
combination  y=0.6  and  C]r=2.5.  It  can  be  observed  that  for  P/8=  1 .0,  Equation 
39  yields  a  value  of  6^=0.091,  hence  in  very  good  agreement  with  proposed 
values  for  this  constant  in  flow  regions  under  local  turbulence  equilibrium 
conditions. 
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Figure  3.4  Variation  of  Cju  with  the  ratio  P/s 
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Limitations  of  Turbulence  Models  Based 
on  Flux-Gradient  Approximations 

As  it  can  be  observed  from  the  expressions  given  in  the  preceding  section, 
most  of  the  assumptions  made  involve  the  use  of  flux  gradient  models.  Over  the 
years,  there  has  been  some  criticism  concerning  the  limitations  of  these  models. 
In  particular,  Corrsin  (1974)  enumerates  some  of  the  necessary  (but  not 
sufficient)  conditions  for  these  assumptions  to  represent  the  actual  processes  in 
terms  of  homogeneity  and  stationarity  of  the  mean  field  being  transported  and  of 
the  turbulence  properties: 

a.  The  transport  mechanism  length-scale  must  be  much  smaller  than  the 
distance  over  which  the  curvature  of  the  mean  transported  field  gradient 
changes  appreciably. 

b.  The  transport  mechanism  time  scale  must  be  much  smaller  than  the  time 
during  which  the  mean  transported  field  changes  appreciably. 

c.  The  transport  mechanism  length  scale  must  be  essentially  constant  over  a 

distance  for  which  the  mean  transported  field  changes  appreciably. 

d.  The  transport  mechanism  velocity  must  be  appreciably  more  uniform  than 
the  length  scale. 

Defining  the  Lagr^uigian  length  scale  for  momentum  transfer,  Z^,  as  the 
product  of  the  Eulerian  (integral)  time  scale  of  the  turbulent  shear  stress,  7^,  and 
the  root-mean-square  of  the  bed-normal  velocity  fluctuations,  w^ms,  whereas  the 
latter  also  is  used  as  a  velocity  scale,  Corrsin  (1974)  expressed  the  former 
requirements  mathematically  as: 


a. 


b.-U 


'%'s 

zt  I  Ts 


l-^l  —  1 
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(40) 

(41) 


c.andd.  U  l^r' 


T  av 


1  dWrms 


I  «  1 


(42) 


where  symbols  like  U^t  represent  second  derivatives  of  mean  velocities  with 
respect  to  the  vertical  coordinate  and  time.  Evaluating  the  former  expressions 
using  boundary  layer  data  from  Blackwelder  and  Kovasznay  (1972),  Corrsin 
found  that  conditions  a  and  b  were  satisfied,  whereas  condition  c  and  d  were 
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violated,  namely  those  requiring  cross-stream  uniformity  of  the  length  scale  and 
root-mean-square  velocity  fluctuations. 

Despite  the  mentioned  violations,  gradient  flux  models  have  surprisingly 
yielded  very  good  agreement  with  experimental  observations  for  a  wide  variety 
of  applications.  As  it  will  be  shown  later,  results  reported  herein  confirm  this 
assertion. 

Numerical  Algorithm 

This  section  will  discuss  the  algorithm  used  for  the  numerical  study  of  steady, 
uniform  open-channel  flow  through  vegetation.  Under  these  conditions  all  the 
differential  equations  used  can  formally  be  reduced  to: 


(43) 


where 

tp  =  any  dependent  variable 

If)  =  associated  exchange  coefficient  defined  as  Tip  — 

Hgff  =  effective  dynamic  viscosity 

=  effective  Prandtl/Schmidt  number 
=  souce  or  sink  term 

Thus  by  assuming  spatial  variations  only  in  the  vertical  direction,  the  numerical 
solution  of  Equation  43  requires  the  discrete  specification  of  V'  in  the  {z,t)  space, 
and  thus  integration  over  the  control  volume  as  shown  in  Figure  3.5.  The  control 
volume  method  proposed  by  Patankar  and  Spalding  (1970)  will  be  used  with  an 
equation  solver  developed  by  Svensson  (1986)  called  PROBE.  A  brief 
description  of  the  numerical  algorithm  follows,  but  the  reader  interested  in  more 
details  is  referred  to  the  aforementioned  references. 
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Figure  3.5  Definition  of  control  volume 
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where  t*  is  some  time  between  U  and  B,  usually  set  equal  to  B  due  to  numerical 
stability  reasons.  Now,  further  decomposition  of  the  terms  in  brackets  in 
Equation  45  in  finite  difference  form  yields: 
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For  simplicity  Equation  46  can  be  rewritten  as: 
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(47) 


zO+j)  B 


«  lT+[ips(i  +  1)  -  ^b(0]  -  7’-[V'b(0  -  V'b(*  -  1)]) 


z(i-h  ^ 


where  =  r^,/Jz(i  +  j)and  T.  =  ./JzO' -  1). 


The  source  term  may  in  turn  be  integrated  as: 


zi‘+j)  B 

// 


Syj  dt  dz 


»  Azil)  S^  f*  At 


(48) 


Furthermore,  it  is  common  to  subdivide  the  source  term  into  two  parts,  one 
containing  the  variable  itself,  as: 

5^  =  5(0  +  5'(0  V'b 


(49) 


so  that  Equation  48  becomes: 


zO'+j)  B 

II 

z(i-h  ^ 
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(50) 


Now,  combining  Equations  44,  47  and  50  yields: 

Az{i)  -  V'£/(0]  «  At  {TA-^PbH  +  1)  -  V'bCO]  -  TAn>B(f)  -  'PbH  -  1)]} 


+  Azil)  At  [5(0  +  5'(0  y>B] 


(51) 


which  may  be  rearranged  as: 


A(0  tpeii  -  1)  +  5(0  V'bCO  +  C(i)  +  1)  =  £>(0 


(52) 
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where: 

AO)  =  T- 


B(i)  =  -  Aii)  -  B(i)  -  ^  +  Az(i)  S\i) 
C{i)  =  T+ 

D(i)  =  -  V't/O  -  ^z(i)  S(i) 


(53) 


It  can  be  observed  that  Equation  52  is  the  compact  expression  for  a 
tri-diagonal  matrix,  and  thus,  once  the  boundary  conditions  are  prescribed. 
Equation  52  can  be  solved  using  for  example  the  Thomas  algorithm  (Patankar 
and  Spalding  1970). 


Regarding  the  boundary  conditions,  basically  two  different  cases  may  be 
distinguished:  (a)  the  value  of  ip  is  prescribed,  or  (b)  the  flux  of  ip  is  given  at  the 
boundary.  In  case  (a)  consider  only  ^(1)  =  ip^  and  ip(N)  =  ipy^,  where  ip^ 
and  Ip  yg  are  the  prescribed  values  of  the  variable  at  the  lower  and  upper 
boundary,  respectively.  For  case  (b)  with  and  y^/  yg  representing  the 
prescribed  fluxes  at  the  lower  and  upper  boundary,  respectively. 


y^LB 


yip,uB 


Ml  +  i) 

f 

Az{N  -  i) 


[ipyiD  -  iPy{\)) 
{iPy(N  -  1)  -  iPyW) 


SO  that: 


(54) 

(55) 


Ml  +  |) 

V'zjd)  =  -Vf.LB  -p - -  +  V'Z)(2) 

Az(.N  -  4) 

V'dW  =  -  y^,c/B  -J - -  +  IpolN  -  1) 


(56) 

(57) 


In  the  particular  case  of  the  k-8  model,  accounting  for  sediment  transport  in 
suspension,  there  will  then  be  a  system  of  four  partial  differential  equations  in 
the  variables  U,  k,  e  and  C,  representing  the  mean  velocity,  turbulent  kinetic 
energy,  rate  of  dissipation,  and  mean  sediment  concentration,  respectively. 
Source  terms  and  boundary  conditions  will  be  treated  as  described  in  the  section 
’’Turbulence  Modeling  by  a  First-Level  Two-Equation  Closure  Scheme”. 
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4  Turbulence  Structure  in 
Open-Channel  Flows 
Without  Vegetation 


Before  the  algorithm  presented  in  the  previous  section  for  simulating  the 
turbulence  stmcture  and  transport  processes  in  vegetated  waterways  is  employed, 
the  capabilities  of  the  model  will  be  tested  for  the  more  often  studied  case  of 
open  channels  without  vegetation.  This  chapter  deals  with  the  verification  of  the 
model  comparing  numerical  predictions  of  the  turbulence  stracture  in 
open-channel  flow,  under  different  roughness  conditions,  against  experimental 
observations  as  well  as  semi-empirical  expressions.  Afterwards,  the  following 
chapter  presents  predictions  of  sediment  transport  processes  in  open  channels 
without  vegetation. 

Mean  Flow 

Mean  velocity  profiles  corresponding  to  two  different  roughness  conditions, 
hydraulically  smooth  and  transitionally  rough  beds,  were  simulated  numerically 
and  the  results  compared  with  the  authors’  experimental  observations.  The 
experiments  were  conducted  under  uniform  flow  conditions  at  the  Hydrosystems 
Laboratory,  University  of  Illinois  at  Urbana-Champaign,  in  a  19.50-m-long, 

0.9 1-m- wide  and  0.61-m-deep  tilting  flume.  Velocity  measurements  were  taken 
with  a  Sontek  acoustic  Doppler  velocimeter  at  a  sampling  frequency  of  25  Hz. 
For  the  smooth-bed  case  the  slope  was  set  to  approximately  0.0006  and  the  mean 
flow  depth  was  0.24  m,  whereas  for  the  transitionally  rough  case  the  slope  was 
about  0.002  with  a  mean  flow  depth  of  0.24  m.  Figure  4.1  compares  model 
predictions  with  the  observations. 
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z 

(m) 


o  ADV  Transitional  (1995) 
-  k-8  Model 


o  ADV  Smooth  (1995) 
-  k-e  Model 


U  (m/s) 


Figure  4.1  Observed  and  predicted  vertical  distributions  of  mean  velocity  for 
experiments  on  transitionally  rough  and  hydraulically  smooth  bed  conditions 

Second-order  Moments 

Dimensionless  values  of  streamwise  and  vertical  standard  deviations  of 
velocity  were  computed,  and  results  were  compared  against  some  of  the  authors’ 
measurements  in  smooth  (Lopez  1994)  and  transitionally  rough  (Nino  1995) 
beds  as  well  as  with  experiments  by  Nezu  (1977)  for  fully  rough  conditions  (the 
subscript  +  indicates  normalization  using  m*  as  scaling  velocity).  Both  the 
authors’  observations  and  Nezu’s  velocity  measurements  were  taken  using 
hot-film  anemometry.  Results  are  shown  in  Figure  4.2  and  Figure  4.3. 

Energy  Budget  Terms 

The  capability  of  the  model  to  simulate  different  terms  in  the  energy  budget 
was  also  checked  by  comparing  experimental  observations  of  turbulent 
production  and  dissipation  rates  with  numerical  results.  Figure  4.4a  illustrates 
the  agreement  for  the  dimensionless  vertical  profile  of  turbulent  production  rate 
in  smooth-bed  flows,  where  data  were  taken  with  the  acoustic  sensor.  Figure 
4.4b  depicts  comparisons  for  turbulence  dissipation  for  smooth  (Lopez  1994)  and 
transitionally  rough  (Nino  1995)  beds,  where  the  solid  line  represents  a 
semiempirical  expression  proposed  by  Nezu  and  Nakagawa  (1993). 

Eddy  Viscosity  and  Mixing  Length 

As  mentioned  in  the  previous  section,  the  boundary  condition  specifying  s  as 
a  function  of  k  at  the  free  surface  allows  for  the  turbulence  damping  at  the 
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a.  Hydraulically  smooth  condition 
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b.  Transitionally  rough  condition 


Figure  4.2  Observed  and  predicted  vertical  distribution  of  dimensionless  rms 
value  of  streamwise  velocity  fluctuations  (Continued) 
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z/H 


c.  Fully  rough  regime 
Figure  4.2  (Concluded) 


Wrms  (m/s) 


Figure  4.3  Observed  and  predicted  vertical  distribution  of  rms  value  of  vertical 
velocity  fluctuations,  smooth-bed  condition 
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a.  Turbulence  production,  smooth-bed  condition 
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b.  Turbulence  dissipation,  smooth-bed  (Lopez  1994)  and  transitionally-rough 
regime  (Nino  1 995) 

Figure  4.4  Observed  and  predicted  vertical  distribution  of  dimensionless 
energy  budget  terms 
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air-water  interface,  i.e.,  reproducing  the  typical  parabolic  profile  of  kinematic 
eddy  viscosity  observed  in  experiments.  Figure  4.5a  and  b  depict  eddy  viscosity 
and  mixing  length  profiles  in  dimensionless  form  for  different  alternatives  in  the 
specification  of  the  boundary  condition,  together  with  experimental  observations 
in  smooth  bed  conditions  conducted  with  the  acoustic  sensor  and  semiempirical 
expressions  proposed  by  Nezu  and  Rodi  (1986). 

Turbulent  Scales 

Dimensionless  vertical  profiles  of  macro  length  scales.  Lx,  were  computed  as: 


L, 


(58) 


where  iTis  a  function  of  Reynolds  number  (Nezu  and  Nakagawa  1993).  Results 
are  illustrated  in  Figure  4.6  together  with  experimental  observations  by  Lopez 
(1994)  and  Nino  (1995). 

Likewise  dimensionless  profiles  of  Taylor  and  Kolmogorov  micro-length 
scales,  X  and  respectively,  were  estimated  as: 


X 

and 


V 


(59) 


(60) 


Figure  4.7a  and  b  compare  the  numerical  results  for  Taylor’s  micro-scale  with 
experimental  observations  and  semiempirical  expressions  by  Nezu  and 
Nakagawa  (1993),  whereas  Figure  4.8a  and  b  show  similar  values  for  the 
Kolmogorov  micro-scale  (therein  Re*  stands  for  the  flow  shear  Reynolds 
number  defined  as  Re*  =  m*  H/v). 

Validity  of  Flux  Gradient  Assumptions 

Following  Corrsin  (1974),  data  collected  with  the  acoustic  sensor  were  used 
to  check  the  validity  of  the  underlying  assumptions  involving  the  use  of  flux 
gradient  models  in  wall-bounded  shear  flows.  The  criteria  represented  by 
Equations  40  and  42  have  been  evaluated  in  Figure  4.9.  As  was  also  observed 
by  Corrsin  (1974)  for  the  case  of  turbulent  boundary  layers,  it  was  found  that  of 


Chapter  4  Turbulence  Structure  in  Open-Channel  Flows  without  Vegetation 


39 


a.  Kinematic  eddy  viscosity 


zlH 


b.  Mixing  length 


Figure  4.5  Observed  and  predicted  vertical  distribution  of  dimensionless 
kinematic  eddy  viscosity  and  mixing  length,  smooth-bed  condition 
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1.0 


Figure  4.6  Observed  and  predicted  vertical  distribution  of  dimensionless 
streamwise  macro-length  scale,  smooth  and  transitionally  rough  bed 
conditions 

all  the  homogeneity  and  stationary  conditions  required  for  the  applicability  of 
gradient  transport  models  in  turbulence,  the  one  requiring  cross-stream 
uniformity  of  the  length  scale  and  root-mean-square  velocity  fluctuations,  i.e. 
Equation  42,  is  the  most  seriously  violated.  However,  as  reflected  by  the  good 
agreement  between  model  predictions  and  experimental  data  in  Figures  4.1  to 
4.8,  the  violation  of  such  requirement  does  not  seem  to  significantly  affect  the 
turbulence  simulation  with  the  model,  at  least  for  the  mean  flow  and  relatively 
low-order  turbulence  statistics  commonly  used  in  engineering  research  (i.e. 
univariate  and  joint  second-order  moments  of  velocity  fluctuations). 
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Suspended  Sediment 
Transport  in 
Open-Channei  Fiows 
Without  Vegetation 


The  capability  of  the  model  for  simulating  (non-cohesive)  suspended 
sediment  transport  in  open  channels  will  be  tested  herein.  The  equation  for  the 
vertical  diffusion  of  sediment  is  solved  together  with  the  momentum,  k-  and 
£-transport  equations,  forming  thus  a  system  of  four  nonlinear  partial  differential 
equations.  This  system  of  equations  becomes  eventually  coupled  if  the 
influence  of  the  sediment  mixture  on  the  flow  structure  is  accounted  for  by 
buoyancy  terms  in  the  equations  for  turbulent  kinetic  energy  and  dissipation 
rates,  or  if  the  density  of  the  two-phase  mixture  becomes  significantly  greater 
than  the  density  of  clear  water. 

Equation  for  Vertical  Sediment  Diffusion 

Vertical  profiles  of  suspended  sediment  concentrations  were  computed  by 
solving  the  equation  for  the  vertical  diffusion  of  sediment,  which  for  uniform 
flow  conditions  reads; 


3  <  C  > 
dt 


w,  <  C  >  -  <  C  w'  >) 


(61) 


where 

<  C  > ,  C  =  spatial/temporal  mean  and  temporal  fluctuating 
suspended  sediment  concentrations,  respectively 

Wj  =  terminal  fall  velocity  of  sediment  particle 

w’  =  bed-normal  temporal  velocity  fluctuation 
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The  first  term  on  the  right  was  approximated  again  using  a  gradient  flux  model, 
<  C  w'  >  =  -  VrlOc  d  <  C  >  jdz,  thus  yielding: 


d  <  C  > 
dt 


Yt  d  <  C  > 
dz  dz 


+  w  ,,  <  C  > 


(62) 


where  Oc  is  the  Prandtl-Schmidt  number  for  sediment  particles.  As  it  will  be 
shown  later,  depending  on  the  assumptions  made,  this  equation  can  be  solved 
either  coupled  or  uncoupled  with  the  set  of  partial  differential  equations  defining 
the  k-s  model.  For  a  given  sediment  particle  size  the  terminal  fall  velocity  was 
estimated  as: 


^v 


4 

3 


(63) 


where 

R  =  submerged  specific  gravity  of  sediment 
Ds  =  mean  sediment  diameter 

Cz).  =  drag  coefficient  of  sediment  particles 

The  drag  coefficient  for  the  sediment  particle  was  computed  with  a  relation  for 
spheres: 


M  [i  +  0.152  +  0.0151  (64) 

with  Rp  =  W.V  D,/v. 

Buoyancy  Effects  upon  Suspended  Sediment 
Transport  Capacity 

The  effect  of  suspended  sediment  (and  its  vertically  variable  concentration 
profile)  upon  the  transport  properties  of  a  stream  is  evidenced  in  three  different 
ways: 

a.  By  affecting  the  turbulence  intensity  of  the  carrier  flow  due  to  the  energy 
spent  in  keeping  the  sediment  suspended.  This  effect  is  commonly 
accounted  for  by  adding  a  buoyancy-related  sink  term  in  the  turbulent 
kinetic  energy  budget  and  a  source  term  in  the  equation  for  the 
dissipation  rate  (Rodi  1984).  Barenblatt  (1953, 1979)  demonstrated  that 
(for  relatively  low  concentrations)  under  similar  flow  conditions  (i.e.  the 
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same  friction  velocity)  a  water  flow  carrying  sediment  in  suspension 
accelerates  under  the  action  of  particles  in  comparison  with  the  clear 
fluid  flow.  This  action  is  mainly  explained  by  a  decrease  in  turbulent 
kinetic  energy,  which  ultimately  causes  a  drag  reduction  (Barenblatt 
and  Golitsyn  1974).  More  recently,  Garcia  (1992)  showed  how 
buoyancy  may  also  be  taken  into  consideration  in  the  boundary 
conditions  for  two-equation  models, as: 


(65) 

k„  =  ^  iPiJiK,,) 

yjCfl 

(66) 

(67) 

where 

^  _  gRw^C„z„ 

=  constant  with  a  value  close  to  5 

and,  as  before,  the  subindex  o  stands  for  the  value  of  the  variable  at  the 
first  grid  point  away  from  the  bed,  which  also  has  to  be  located  within 
the  equilibrium  layer.  It  can  be  shown  (Lopez  1997)  that  the  former  two 
approaches  yield  identical  results,  with: 


^  (68) 

(69) 

(70) 

and 

^  ^  X  K*/w,-l  (71) 


It  would  be  of  interest  to  test  if  the  value  of  as  prescribed  by  equation 
71,  is  indeed  constant  as  suggested  by  observations  of  atmospheric 
boundary-layer  flows. 

b.  By  changing  the  density  of  the  mixture,  Qm'- 

(1  -h  /?  Q  (72) 

with  Qw  representing  the  density  of  clear  water.  Or  by  changing  the 
viscosity  of  the  mixture,  v^,  where  according  to  Einstein  (Graf  1971): 

=  v„il+keC)  (73) 


Chapter  5  Suspended  Sediment  Transport  in  Open-Channel  Flows  without  Vegetation 


47 


with  denoting  the  viscosity  of  clear  water.  Happel  and  Brenner 
(1965)  found  the  Einstein  viscosity  constant,  ke,  to  be  =  2.5. 


c.  By  changing  the  Prandtl-Schmidt  number,  Oc  (Launder,  Reece  and  Rodi 
1975),  in  the  sediment  diffusion  equation. 


1  +  0,'(C/  -  <P,)  B 
I  +4>  0,B 


(74) 


where 

(Jco  =  the  Prandtl/Schmidt  number  under  nonstratified 

conditions 

(pc,  (pc',  Cc  =  model  parameters,  that  for  the  case  of 

temperature-induced  stratification  are  equal  to 
0.31, 0.16  and  1.60,  respectively 

B  =  dimensionless  buoyancy  parameter  defined  as  _  or!^^ 

*  dz 

The  influence  of  each  of  these  factors  depends  upon  the  particular  problem 
under  consideration.  Several  experimental  (Vanoni  1946;  Einstein  and  Chien 
1955;  Coleman  1981)  and  theoretical  (Barenblatt  and  Golitsyn  1974)  works  have 
shown  how  the  velocity  of  the  flow  increases  with  the  mean  sediment 
concentration.  However  the  interaction  of  suspended  matter  and  the  turbulence 
structure  of  the  flow  is  to  date  not  fully  understood.  In  the  present  work  seven 
different  combinations  of  the  factors  discussed  in  the  preceding  paragraph  were 
numerically  analyzed  in  order  to  determine  the  effects  on  the  suspended  sediment 
transport  capacity  of  the  flow  and  finally  decide  which  model  to  use  in  the 
presence  of  vegetation.  Combination  1  corresponded  to  a  decoupled  solution  of 
the  vertical  diffusion  equation,  therefore  neglecting  all  the  factors  mentioned. 
Combination  2  dealt  with  the  inclusion  of  the  buoyancy-related  terms  in  the 
transport  equations  for  k  and  s,  but  maintaining  the  standard  wall  functions  and 
keeping  constant  both  the  density  of  the  mixture  and  Oc-  In  combination  3, 
buoyancy  was  considered  and  the  density  of  the  mixture  was  allowed  to  vary 
vertically  with  the  concentration.  Combination  4  corresponded  to  the  inclusion 
of  buoyant  terms  plus  variations  in  density  and  Oc.  In  the  fifth  combination, 
density  and  ot  were  kept  constant  while  considering  both  buoyant-related  terms 
and  modified  wall  functions.  Combination  6  is  similar  to  the  former  except  that 
density  was  also  allow  to  vary  with  concentration.  And  lastly,  the  seventh 
combination  considered  both  variations  in  density  of  the  mixture  and  Oc  while 
also  introducing  the  buoyancy-related  terms  and  the  modified  wall  functions. 


Basically  two  different  outputs  were  considered,  vertical  profiles  of 
suspended  sediment  concentration  and  total  transport  capacity,  obtained  by 
vertically  integrating  the  product  of  mean  velocity  and  local  sediment 
concentration.  Results  of  vertical  concentration  profiles  were  compared  against 
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the  Rousean  distribution: 


_C 

C, 


(H-z)  jH-b) 
z  b 


HU, 


(75) 


where 

b  =  0.05  H 

Q,  =  near-bed  sediment  concentration 

The  sediment  concentration  near  the  bed  was  estimated  using  the  expression  of 
Garcia  and  Parker  (1991): 


AZl 

o3o  ^“) 


(76) 


with  Z„  =  u^/w,  Rep^-^  2a\d  A  =  1.30  .c  10  Equation  76  provides  a  mean 

of  estimating  the  bed  sediment  concentration  under  equilibrium  conditions  at  z/H 
=  0.05,  and  was  derived  using  data  covering  the  following  ranges  of  the 
variables: 


2.0  10"^  <  Cf,<  6.0  10"2; 

0.70  <  u^/ws,  <  7.50; 

240  <  H/D,  <  2,400; 

3.50  <  Rep  <  37.0. 

The  computed  suspended  transport  capacity,  q^s,  was  compared  against 
predictions  from  the  formula  due  to  Einstein  (1950): 


11.57  Cj.  M.  b 


(77) 


where  kc  represents  a  measure  of  the  roughness,  and  /j  and  I2  are  numerically 
evaluated  integrals  (see  also  Graf  1971). 


It  must  be  said  that  a  particular  limitation  arises  when  trying  to  solve  all  the 
equations  as  a  coupled  system.  This  limitation  is  based  on  the  fact  that  the  wall 
functions  have  to  be  evaluated  at  a  grid  point,  Zo,  located  between  30  and  100 
wall  units  from  the  bed,  whereas  the  proposed  expression  for  estimating  the 
near-bed  concentration  gives  the  value  of  Q,  at  2;,  =0.0577.  Therefore,  for  a 
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smooth-bed  two-dimensional  open  channel: 


0.05H  Jg  So  H 
30  <  - —  <  100 

where  So  =  bed  slope.  Equation  78  constitutes  an  important  constraint  that 
specifies  the  required  range  of  H  for  a  given  slope  and  water  temperature. 


(78) 


Suspended  sediment  transport  capacity  was  computed  for  a  two-dimensional 
open  channel  with  So  =  0.001  and  H  =  0.07  m  for  different  sediment  sizes  using 
the  seven  aforementioned  combinations.  Results  are  depicted  in  Figure  5.1  in 
dimensionless  form  together  with  predictions  from  the  Einstein’s  (1950)  formula. 
In  order  to  appreciate  the  effects  of  buoyancy  on  water  discharge.  Figure  5.2 
illustrates  the  ratio  qw^qwo  corresponding  to  the  same  conditions  as  Figure  5.1, 
where  is  the  computed  water  discharge  for  each  combination  and  q^o  is  the 
computed  water  discharge  for  combination  1,  thus  without  considering  any 
sediment-turbulence  interaction. 

These  graphs  show  that,  while  water  discharge  increases  as  much  as  56  percent 
depending  upon  the  sediment  concentration  and  calculation  procedure, 
suspended  sediment  transport  capacity  is  relatively  well  predicted  by  neglecting 
sediment-turbulence  interactions,  especially  for  low  concentrations.  Based  on 
these  results  and  the  fact  that  the  sediment-laden  flows  considered  herein  have 
relatively  low  sediment  concentrations,  the  results  shown  in  the  next  two 
paragraphs  were  computed  neglecting  buoyancy  effects. 

Vertical  Profile  of  Suspended  Sediment 
Concentration 


Figure  5.3  shows  computed  dimensionless  vertical  profile  of  suspended 
sediment  concentration  together  with  predictions  by  the  Rousean  model  for  three 
different  sediment  sizes,  namely  40, 100  and  150  pm  (Rep  =  1.87, 4.02  and 
11.38,  respectively). 

Suspended  Sediment  Transport  Capacity 

Figure  5.4  shows  results  of  the  variation  in  suspended  sediment  transport 
capacity,  for  a  two-dimensional  channel  {Sg  =  0.0036  and  H  =  0.35  m)  for 
seven  different  sediment  sizes,  whereas  Figure  5.5  depicts  this  capacity  as  a 
function  of  flow  depth  for  given  values  of  Sg  and  Dg.  Results  of  the  numerical 
model  are  compared  in  both  figures  to  predictions  with  the  formula  by  Einstein 
(1950). 
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Figure  5.1  Suspended  sediment  transport  capacity  as  function  of  Rouse 
number  and  different  buoyancy  effects  together  with  predictions  by  Einstein’s 
(1 950)  formula 
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Figure  5.2  Relative  water  discharge  as  function  of  Rouse  number  and 
different  buoyancy  effects 
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Figure  5.3  Dimensionless  vertical  profile  of  suspended  sediment  concentration 
for  different  mean  diameters 
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Figure  5.4  Estimated  suspended  sediment  transport  capacity  for  different 
mean  diameters 
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Figure  5.5  Estimated  suspended  sediment  transport  capacity  for  different 
mean  flow  depths 
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6  Turbulence  Structure  in 
Vegetated  Open-Channel 
Flows 


In  this  chapter  experimental  results  for  open-channel  flows  in  the  presence  of 
both  rigid  and  flexible  simulated  vegetation  will  be  used  to  check  the  validity  of 
the  assumptions  made  previously.  First,  the  experimental  conditions  will  be 
summarized.  Then,  computed  vertical  profiles  of  mean  flow  as  well  as 
turbulence  characteristics  will  be  compared  against  experimental  results.  Lastly, 
the  impact  of  using  wall  functions  different  from  the  usual  ones  in  the  standard 
model  will  be  briefly  explored  and  Corrsin’s  criteria  for  the  applicability  of 
gradient-flux  models  in  the  presence  of  vegetation  will  be  evaluated. 

Experiments  in  Open-Channel  Flows 
with  Simulated  Vegetation 

As  mentioned  in  Chapter  3,  a  series  of  experiments  was  conducted  at  the 
Hydrosystems  Laboratory  with  the  goal  of  characterizing  the  drag  coefficient  in 
the  particular  case  of  free-surface,  turbulent  flows  and  of  providing  information 
for  the  verification  of  the  k-e  model.  Regarding  the  second  reason,  it  was 
considered  crucial  for  the  authors  to  conduct  their  own  observations  because 
most  of  the  data  available  lacked  a  detailed  description  of  the  experiments;  in 
particular,  no  clear  specification  of  the  measuring  location  and  the  spatial 
averaging  procedure  are  given.  From  the  discussions  in  Chapter  3,  the 
importance  of  the  averaging  procedure  employed  to  determine  one-dimensional 
parameters  becomes  clear. 

The  experiments  were  conducted  under  uniform  flow  conditions  in  a 
19.50-m-long,  0.91-m-wide  and  0.61-m-deep  tilting  flume  (Dunn,  Ldpez  and 
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Flexible 


Garcia  1996).  Velocity  measurements  were  taken  with  a  Sontek  acoustic 
Doppler  velocimeter  at  a  sampling  frequency  of  25  Hz,  at  four  different 
plan-locations  with  10  measuring  points  in  each  vertical.  Cylindrical  wooden 
dowels  and  commercial  drinking  straws  were  used  to  simulate  rigid  and  flexible 
vegetation,  respectively.  Table  6.1  provides  a  summary  of  the  experimental 
conditions,  where  a  represents  plant  density,  Sq  is  the  averaged  bed 
(water-surface)  slope,  Q  is  the  discharge,  H  represents  the  normal  flow  depth,  a 
is  the  averaged  deflection  angle  of  the  simulated  plants.  Re  is  the  Reynolds 
number  and  Fr  is  the  Froude  number. 

TABLE  6.1  Experimental  Conditions 


Exp# 

a 

.^0 

Q 

H 

a 

Re 

Fr 

(l/m) 

(m3/s) 

(m) 

n 

Exp.  1 

1.09 

0.0036 

0.179 

0.335 

0.00 

224,000 

0.33 

Exp.  2 

1.09 

0.0036 

0.088 

0.229 

0.00 

113,000 

0.29 

Exp.  3 

1.09 

0.0036 

0.046 

0.164 

0.00 

57,000 

0.24 

Exp.  4 

1.09 

0.0076 

0.178 

0.276 

0.00 

191,000 

0.36 

Exp.  5 

1.09 

0.0076 

0.098 

0.203 

0.00 

125,000 

0.37 

Exp.  6 

0.27 

0.0036 

0.178 

0.267 

0.00 

196,000 

0.39 

Exp.  7 

0.27 

0.0036 

0.095 

0.183 

0.00 

120,000 

0.42 

Exp.  8 

2.46 

0.0036 

0.180 

0.391 

0.00 

258,000 

0.29 

Exp.  9 

2.46 

0.0036 

0.058 

0.214 

0.00 

69,700 

0.19 

Exp.  10 

2.46 

0.0161 

0.180 

0.265 

0.00 

203,000 

0.40 

Exp.  11 

0.62 

0.0036 

0.177 

0.311 

0.00 

222,000 

0.35 

Exp.  12 

0.62 

0.0110 

0.181 

0.233 

0.00 

238,000 

0.58 

Exp.  13 

1.09 

0.0036 

0.179 

0.368 

35.0 

228,000 

0.28 

Exp.  14 

1.09 

0.0101 

0.180 

0.232 

51.0 

257,000 

0.62 

Exp.  15 

1.09 

0.0036 

0.093 

0.257 

34.0 

112,000 

0.23 

Exp.  16 

0.27 

0.0036 

0.179 

0.230 

65.0 

227,000 

0.56 

Exp.  17 

2.46 

0.0036 

0.078 

0.279 

12.0 

94,900 

0.18 

Exp.  18 

2.46 

0.0101 

0.179 

0.284 

45.0 

250,000 

0.45 

Mean  Flow  Characteristics 

Model  predictions  regarding  mean  (spatial  and  temporal  averaged)  velocities 
were  compared  against  the  experimental  observations.  The  modeled  turbulent 
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kinetic  energy  in  the  two-equation  turbulence  model  corresponded  to  Equation 
10;  hence  values  of  =  1.0  and  =  1.33  were  used  in  the  computations. 
Figure  6.1  compares  numerical  simulations  and  experimental  observations  for 
two  of  the  rigid- vegetation  tests  in  Table  6.1,  whereas  Figure  6.2  shows 
predictions  corresponding  to  the  flexible  conditions.  As  will  be  also  observed  in 
the  Reynolds  stress  computations,  the  existence  of  secondary  currents  seems  to 
play  an  important  role  above  the  simulated  vegetation,  both  in  retarding  the  flow 
and  decreasing  the  turbulent  momentum  transfer  in  the  vertical.  Hence,  since 
these  computations  simulate  two-dimensional  flows,  computed  mean  velocities 
above  the  plant  canopy  are  slightly  larger  than  the  measured  ones. 

Figure  6.3  depicts  the  vertical  mean  velocity  profile  of  two  experiments  (rigid 
and  flexible  conditions)  in  semilog  scale,  showing  the  agreement  between  both 
numerical  and  experimental  results  with  the  logarithmic  law.  The  differences  in 
slope  between  the  model  predictions  and  the  measured  profile  may  be  explained 
by  the  different  values  of  the  predicted  and  observed  turbulent  momentum 
transfer  close  to  the  top  of  the  simulated  canopy  (see  ’’Reynolds  stresses”  in  next 
section).  Indeed,  the  vertical  slope  of  the  velocity  profile  in  the  equilibrium 
layer  of  wall-bounded  flows  becomes  dU jdz  —  z,  which  clearly  explains 
that  a  larger  value  of  the  shear  velocity  yields  larger  slopes  at  the  same  distance 
from  the  bed.  Note  also  how  the  numerical  model  predicts  the  existence  of  a 
region  immediately  above  the  simulated  plants,  where  the  velocity  is  not 
logarithmically  distributed,  in  agreement  with  reported  results  on  velocity 
distributions  in  the  roughness  sublayer  (Lopez  1997). 

Second-order  Moments 

This  section  shows  the  capability  of  the  numerical  model  to  simulate  the 
vertical  stmcture  of  second-order  moments,  i.e.  Reynolds  stresses  and 
turbulence  intensities. 

Reynolds  stresses 

Computed  vertical  profiles  of  spatially  averaged  Reynolds  stresses  are 
depicted  in  Figure  6.4  and  6.5  for  rigid  and  flexible  conditions,  respectively, 
together  with  experimental  observations,  for  the  same  experiments  as  Figure  6.1. 
As  mentioned  previously,  a  very  good  agreement  between  experimental  values 
and  model  predictions  is  observed  for  flow  within  the  simulated  vegetation, 
whereas  the  measured  Reynolds  stresses  above  the  simulated  canopy  were 
consistently  smaller  than  the  theoretical  ones  for  two-dimensional  open-channel 
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Figure  6.1  Observed  and  predicted  vertical  distribution  of  mean  velocity,  rigid 
conditions 
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Figure  6.2  Observed  and  predicted  vertical  distribution  of  mean  velocity, 
flexible  conditions 
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Expl 


a.  Rigid  Condition  U  (mis) 


U  (m/s) 

b.  Flexible  condition 

Figure  6.3  Observed  and  predicted  vertical  distribution  of  mean  velocity  in 
semilog  scale 
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Figure  6.4  Observed  and  predicted  vertical  distribution  of  Reynolds  stresses, 
rigid  conditions 


Chapter  6  Turbulence  Structure  in  Vegetated  Open-Channel  Flows 


-  <  u'w'  >  (m^ls^) 


Figure  6.5  Observed  and  predicted  vertical  distribution  of  Reynolds  stresses 
flexible  conditions 
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flows.  However,  this  phenomenon  is  typical  for  free-surface  flows  (Nezu  and 
Nakagawa  1993),  the  deviation  being  explained  by  the  action  of  secondary 
currents  as  well  as  other  components  of  the  Reynolds  stress  tensor,  the 
magnitude  of  this  effect  being  a  function  of  the  width-to-depth  ratio  (aspect 
ratio). 

Turbulent  kinetic  energy  and  turbulence  intensities 

Computations  of  experimental  values  of  total  turbulent  kinetic  energy  and 
turbulence  intensities,  <  m/'  m,"  > ,  are  very  difficult  to  obtain,  especially  due 
to  the  large  number  of  measuring  locations  needed  to  obtain  representative 

values  of  <  uj  uj  > .  It  is  worth  noting  again  that  for  i  ^  j  (i.e.  off-diagonal 
components  of  the  total  turbulent  stress  tensor)  the  contribution  of  the 
wake-related  production  term  is  negligible  (Raupach  et  al.  1986).  Hence,  in 
order  to  compare  numerical  results  with  experimental  observations  of 
<  m/  m,'  > ,  the  model  was  run  with  =  0.0  and  =  0.0.  Figures  6.6 
and  6.7  depict  the  results  obtained  for  the  streamwise  turbulence  intensities 
together  with  the  experimental  data.  In  order  to  better  visualize  the  difference 
between  <  m/  h/  >  and  <  u/'  u/'  >,  Figure  6.8  illustrates  the  computed 
values  of  the  total  streamwise  intensities  obtained  when  the  model  was  ran  using 
Cjj.  =  l.OandC^g  =  1.33. 

In  light  of  these  results  and  the  discussion  in  chapter  3,  it  becomes  clear  why 
very  small  weighting  coefficients  in  the  drag-related  source  terms  for  the  k  and  e 
equations  yield  very  good  predictions  of  the  observed  values  of  turbulence 
intensities  in  water  flows. 

Energy  Budget  Terms 

To  further  investigate  the  performance  of  the  model,  different  terms  in  the 
turbulent  kinetic  energy  budget  were  computed,  and  when  possible  compared 
against  the  experimental  observations. 

Spatially  averaged  time-mean  values 

As  mentioned  in  the  previous  section,  only  accurate  measurements  of 
spatially  averaged  time-mean  values  could  be  obtained  from  the  experiments, 
and  those  are  the  results  that  are  simulated  herein.  Figure  6.9  and  6.10  illustrate 
vertical  profiles  of  the  different  terms  in  the  energy  budget  made  dimensionless 
using  U*hp  and  hp  as  scaling  velocity  and  length  scales,  respectively.  U*hp  is  the 
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Figure  6.7  Observed  and  predicted  vertical  distribution  of  streamwise 
turbulence  intensity,  flexible  condition 
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a.  Rigid  vegetation 


y<  >  (mis) 


b.  Flexible  vegetation 

Figure  6.8  Observed  and  predicted  vertical  distribution  of  total  streamwise 
turbulence  intensity 
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dt  ^*hc 
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dt  ^*hc 

Figure  6.9  Observed  and  predicted  vertical  distribution  of  different  terms  in  the 
spatially  averaged,  temporal-mean,  turbulent  kinetic  energy  budget  for 
(Cfi(,Cfe)=(0.0,0.0).  Lines  represent  model  predictions  and  symbols  are 
observed  values. 
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Figure  6.10  Observed  and  predicted  vertical  distribution  of  different  terms  in 
the  spatially  averaged,  temporal-mean,  turbulent  kinetic  energy  budget  for 
(Cfi(,Cfe)=(0.25,0.33).  Lines  represent  model  predictions  and  symbols  are 
observed  values. 
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square  root  of  the  Reynolds  stress  at  the  top  of  the  simulated  canopy  and  hp 
represents  the  average  plant  height. 

Budget  of  total  turbulent  kinetic  energy 

Although  no  experimental  observations  were  available  for  comparison, 
dimensionless  vertical  profiles  of  the  total  turbulent  kinetic  energy  were 
computed  (i.e.,  =  1.0  and  =  1.33).  Figure  6.11  depicts  the  results 

obtained  for  the  same  experiments  as  in  Figure  6.1. 

Eddy  Viscosity  and  Mixing  Length 

Again,  only  reliable  values  of  spatially  averaged  time-mean  values  of  eddy 
viscosity  and  mixing  length  could  be  obtained  from  the  experimental 
measurements.  Figures  6.12  to  6.14  compare  experimental  observations  with 
numerical  results  using  three  different  sets  of  values  for  the  weighting  factors  of 
the  drag-related  terms,  namely  (C^,  Cj^  =  (0.0,  0.0), 

(C^,  Cp  =  (0.8,  1.04),  and  (C^,  Cp  =  (1.00,  1.33). 

Turbulent  Length  Scales 

Dimensionless  vertical  profiles  of  macro-  and  micro-length  scales  were 
computed  using  Equations  58,  59  and  60.  Numerical  results  corresponding  to 
the  conditions  of  Exp  1  are  plotted  against  experimental  observations  in  Figure 
6.15. 

Momentum  Transfer  to  the  Bed 

In  the  evaluation  of  suspended  sediment  transport  processes  in  open  channels, 
the  accurate  estimation  of  the  momentum  transfer  to  the  bed  plays  a  crucid  role, 
since  this  value  is  typically  used  to  evaluate  the  ability  of  the  flow  to  entrain 
sediment  from  the  bed.  Figure  6.16  illustrates  the  variation  of  the  shear  velocity, 
defined  as  the  square  root  of  the  bed-shear  stress  per  unit  density,  as  a  function  of 
normalized  plant  density.  Note  that  in  this  graph  the  shear  velocity  has  been 

standarized  using  the  total  streamwise  momentum  due  to  gravity,  i.e.  /gHSo. 

Manning’s  Resistance  Coefficient 

In  order  to  evaluate  the  effect  of  vegetation  upon  flow  resistance,  values  of 
Manning’s  resistance  coefficients  were  computed  as: 
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Z)(<  u-  m/  >)  hp 
dt  ^*hc 


Figure  6.11  Observed  and  predicted  vertical  distribution  of  different  terms  in 
the  total  turbulent  kinetic  energy  budget 
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Figure  6.13  Observed  and  predicted  vertical  distribution  of  mixing  length,  rigid 
condition 
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b.  Mixing  length 

Figure  6.14  Observed  and  predicted  vertical  distribution  of  eddy  viscosity  and 
mixing  length,  flexible  vegetation 
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Figure  6.15  Observed  (symbols)  and  predicted  (thin  lines)  vertical  distribution 
of  dimensionless  length  scales  for  Expl 
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hH  So 


Figure  6.16  Momentum  transfer  towards  the  bed  vs.  dimensionless  plant 
density 

//5/3  s]/^  (79) 

"  “  qw 
where 

qw  =  specific  water  discharge 

Figure  6.17  shows  the  variation  of  n  with  plant  density,  a,  for  given  values  of 
water  discharge,  =  1.12  revivals,  channel  slope,  Sg  -  0.0036,  and  plant  height, 
hp  =  0.10  m.  Since  in  engineering  practice  it  is  more  common  to  measure  density 
as  number  of  plants/stems  per  square  meter,  X^,  Figure  6.18  depicts  the  variation 
of  both  Manning’s  n  and  flow  depth  with  X^  for  D  =  6.4  mm  (diameter  of 
cylinders  used  in  authors’  experiments).  As  can  be  observed,  the  resistance 
coefficient  remains  almost  constant  for  low  densities  and  it  shows  a  sharp 
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Figure  6.17  Computed  values  of  Manning’s  resistance  coefficient  as  a  function 
of  plant  density 
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Figure  6.1 8  Computed  values  of  flow  depth  and  Manning’s  resistance 
coefficient  as  a  function  of  density  measured  as  number  of  plants 
per  square  meter 
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increase  after  a  critical  density  value  has  been  reached,  increasing  linearly 
afterwards.  As  mentioned  in  Chapter  2,  this  linear  increase  has  also  been 
observed  in  the  field  (Freeman,  Hall  and  Abraham  1994). 


Impact  of  Wall  Functions 


As  mentioned  in  Chapter  2,  the  wall  functions  used  as  bottom  boundary 
conditions  in  the  standard  k-e  model  implicitly  assume  the  existence  of  an 
equilibrium  layer  (match  or  overlap  region),  where  the  logarithmic  law  describes 
the  vertical  distribution  of  mean  velocities.  As  can  be  clearly  observed  from  the 
budgets  of  turbulent  kinetic  energy  this  does  not  hold  for  the  one-dimensional 
description  of  flow  through  vegetation.  However,  it  was  observed  that  the  exact 
form  of  the  wall  functions  was  of  little  relevance  for  the  computations,  and  that 
the  drag-related  source  terms  were  of  more  critical  importance  for  the  flow 
structure.  As  an  example.  Figure  6.19  shows  two  different  computations 
corresponding  to  Expl,  where  the  wall  functions  at  the  bed  were  modified  as 
follows: 


U„ 


=  C 


u* 
vel  X 


ko 


=  Ck 
=  Ce 


Zo 


(80) 

(81) 

(82) 


Indeed,  specified  values  of  mean  velocities  and  shear  stress  using  these 
expressions  are  so  small  that  the  simulated  turbulence  stracture  adjusts  itself  to  a 
common  profile,  with  results  almost  insensitive  to  the  exact  value  of  the 
boundary  conditions  at  the  bed.  In  view  of  these  results,  the  standard  wall 
functions  have  been  used  throughout  the  present  work,  i.e.,  =  1.0, 

Cj^  =  1.0  and  Ce  =  1.0. 


Validity  of  Gradient-Flux  Assumptions 
for  Flow  through  Vegetation 

Following  the  same  ideas  as  in  Chapter  4,  data  collected  with  the  acoustic 
sensor  were  used  to  check  the  validity  of  the  underlying  assumptions  involving 
the  use  of  flux  gradient  models  in  vegetated,  free-surface,  wall-bounded  shear 
flows.  The  criteria  represented  by  Equations  40  and  42  have  been  evaluated  in 
Figure  6.20,  where  results  corresponding  to  flow  without  vegetation  are  also 
presented  for  comparison.  Again,  it  may  be  clearly  observed  that  from  all  the 
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•  With  Vegetation 
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required  conditions,  the  one  requiring  cross-stream  uniformity  of  the  length 
scale  and  root-mean-square  velocity  fluctuations,  i.e..  Equation  42,  is  the  most 
seriously  violated.  As  evidenced  in  previous  figures,  however,  the  violation  of 
this  requirement  does  not  seem  to  significantly  affect  the  turbulence  simulations 
of  the  model,  at  least  for  the  statistics  commonly  used  in  engineering  research. 


Final  Remarks 

The  graphs  presented  above  demonstrate  the  overall  ability  of  the  numerical 
model  to  simulate  not  only  the  most  commonly  used  flow  statistics,  like  mean 
velocity,  turbulence  intensities  and  Reynolds  stress,  but  also  different  terms  in 
the  budget  of  turbulent  kinetic  energy  as  well  as  mixing  properties  and 
turbulence  macro  and  micro  length  scales  for  flow  through  vegetation.  As 
mentioned  before,  observed  profiles  of  temporal  mean  variables  averaged  over 
space  are  being  best  represented  using  negligible  values  for  the  drag-related 
weighting  coefficients.  Differences  between  numerical  results  using  values  of 
(C^,  C^)  =  (0.0,  0.0)  and  (Cyj,  C^)  =  (1.00,  1 .33)  are  found  to  be  larger 
at  the  top  of  the  simulated  plants,  and  to  decrease  towards  the  free  surface.  If 
real,  these  differences  would  imply  the  existence  of  streamwise  turbulence 
heterogeneities  in  flow  regions  above  the  plants.  More  research  is  however 
needed  to  particularly  address  the  influence  of  flow  and  vegetation  properties 
upon  the  ratio  between  turbulent  micro-length  scales  and  element  wakes,  i.e.  the 
relationship  between  these  properties  and  the  values  of  the  weighting  coefficients 
in  the  drag-related  terms.  Hereafter  all  the  one-dimensional  (thus  involving  the 
use  of  spatial  and  temporal  averaged  conservation  laws)  numerical  simulations  of 
turbulence  processes  will  be  performed  with  values  for  these  coefficients  of 
C^)  =  (1.00,  1.33),  i.e.  using  Equation  10  for  the  turbulent  kinetic 

energy  balance  and  therefore  k  =  (<  uj'  uj  '  >  +  <  ul  u  ’  >)/2. 
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7  Suspended  Sediment 
Transport  in  Vegetated 
Open  Channels 


In  previous  chapters  two-equation  turbulence  model  predictions  have  been 
verified  against  experimental  observations  both  in  open  channels  and  vegetated 
free-surface  flows,  and  the  capability  of  the  algorithm  for  simulating  the 
turbulence  structure  has  been  checked.  The  last  stage  of  the  investigation  will  be 
then  to  apply  the  numerical  code  to  estimate  the  suspended  sediment  transport 
capacity  of  flows  in  vegetated  waterways.  The  present  section  begins  with  the 
introduction  of  dimensional  analysis,  which  allows  for  the  identification  of  the 
dimensionless  parameters  that  govern  the  problem.  Afterwards,  the  data  set  by 
Tollner,  Barfield  and  Hayes  (1982)  will  be  employed  to  check  the  outcome  of  the 
numerical  experiments.  Further  results  include  comparisons  between  simulated 
vertical  distributions  of  sediment  concentration  and  predicted  profiles  using  the 
Rousean  distribution,  estimations  of  transport  capacity  as  function  of  different 
dimensionless  parameters,  and  computations  of  relative  transport  capacity 
between  vegetated  and  nonvegetated  open  channels  under  similar  hydraulic 
conditions,  as  a  function  of  plant  density,  sediment  diameter,  etc. 


Dimensional  Analysis  of  Sediment  Transport 
in  Vegetated  Open  Channels 

The  investigation  of  sediment  transport  processes  in  vegetated  channels 
involves  the  consideration  of  so  many  variables  characterizing  the  sediment, 
flow  and  plants  properties,  that  the  problem  might  seem  almost  intractable. 
Dimensional  analysis  constitutes  a  valuable  tool  in  these  situations.  Any 
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variable,  x,  characterizing  sediment  transport  processes  in  vegetated  waterways 
may  be  expressed  as  a  function  of  the  same  variables  that  govern  the 
phenomenon  in  open  channels  without  vegetation,  plus  some  new  variables  that 
characterize  the  plant  properties: 

X  =  /[“*’  (Qs-Q),  Ds,  D,  a,  hp,  a]  (83) 

with  a  representing  a  dimensionless  parameter  defining  the  flexibility  of  the 
vegetation.  It  is  worth  noting  that  in  the  previous  expression  m*  is  the  bed-shear 
velocity  associated  with  the  average  streamwise  momentum  transfer  to  the  bed. 


If  M*,  ^  and  are  selected  as  the  basic  quantities  for  performing  the  standard 
dimensional  analysis  (see  for  example  Yalin  1977),  then: 


y  =  R  rR  D  a^a\ 

where  x  is  the  dimensionless  form  of  the  transport  property  and 
R  =  (Qs  ~  Q)/Q-  It  may  be  further  observed  that: 


(84) 
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Dsjg  R  Ds 


h  R  Ds 


Jg  R  R>s 
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g  R  Ds  w 


(85) 


(86) 


g  R  Ds 

Moreover,  it  can  be  shown  that  for  grains  of  a  given  shape,  the  dimensionless 
parameter  w]/{gRD^  can  be  represented  solely  as  a  function  of  Rep  (Parker  1978; 
Dietrich  1982).  With  all  these  considerations  plus  some  additional  combinations 
of  dimensionless  parameters.  Equation  84  may  be  rewritten  as: 


X  H^^P'Ws’Ds'^'Ds’^  ^ 


hp  \ 


(87) 


Generally,  calculations  will  be  based  on  constant  values  of  /?  =  1 .65  and  very 
large  values  of  the  ratio  H/Ds  and  D/D^;  therefore  Equation  87  may  be  further 
simplified  into: 


R. 


Z.  T4  n  — il  n 


(88) 


Experiments  by  Tollner,  Barfield  and  Hayes  (1982) 

As  mentioned  in  the  introduction,  there  are  not  many  reliable  data  sets 
available  with  suspended  sediment  information  for  vegetated  open  channels. 

The  complexity  of  the  problem  led  some  investigators  to  study  sediment 
transport  processes  in  the  laboratory  by  simulating  vegetation  with  different 
elements.  In  this  section  the  predictions  of  the  k-e  model  will  be  compared  with 
the  experimental  observations  by  Tollner,  Barfield  and  Hayes  (1982).  These 
investigators  used  a  relatively  narrow  and  short  laboratory  flume  (0.13-m-wide, 
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O.lO-m-deep  and  2.10-m-long)  with  16d  finish  nails  (Tollner,  private 
communication)  simulating  rigid  vegetation  in  a  staggered  pattern.  They 
simulated  sediment  using  glass  beads  of  five  diffferent  mean  diameters.  They 
found  that  their  experimental  results  for  suspended  sediment  transport  could  be 
well  described  using  a  modified  version  of  Graf’s  parameters,  as: 


,  -0.153 

W  =  0.52  (^f) 

where: 


(89) 


=  Co  <  Ml  > 
and 


W  = 


R  Ds 

Rs 


where 


(91) 


Co  =  the  volumetric  suspended  load  concentration 
Rs  =  an  “equivcilent  hydraulic  radius”  defined  as: 

R  =  ^  (92) 

^  bp  +  2  H 

with  bp  representing  the  plant  spacing. 


It  is  worth  mentioning  that  Equation  89  does  not  contain  all  the 
dimensionless  parameters  specified  in  Equation  88,  in  particular  the  parameters 
Ha  and  H/hp  are  missing.  The  absence  of  the  latter  is  justified  because  the 
authors  seemed  to  have  reported  only  values  corresponding  to  emergent 
vegetation.  The  absence  of  the  parameter  Ha  may  also  be  justified  by  the  fact 
that  flow  depth  and  ntdl  spacing  varied  over  a  narrow  range  (the  latter  in  the 
range  0.945-1.583  cm). 


In  order  to  compare  numerical  results  with  experimental  observations  the  four 
hydraulic  conditions  given  in  Tollner  (1974)  were  selected,  where  some 
turbulence  measurements  were  conducted  as  well  using  hot-film  anemometry. 
Figure  7.1  illustrates  the  results  of  the  k-e  model  compared  to  predictions  by 
Equation  89.  As  can  be  observed,  experimental  results  corresponding  to  the 
higher  slope  are  very  well  predicted  by  the  numerical  model,  while  computed 
sediment  transport  capacity  for  the  smallest  slope  is  higher  than  the  observed 
values.  The  observed  disagreement  may  be  attributed  to  the  small  dimensions  of 
the  experimental  facilities,  which  may  have  precluded  the  establishment  of 
equilibrium  conditions.  The  shallow  flow  depths  used  in  the  experiments, 
ranging  from  about  1.3  to  5.0  cm,  suggest  also  that  the  flows  might  not  have  had 
large  enough  Reynolds  numbers,  while  the  numerical  model  developed  in  the 
present  work  applies  only  to  high-Reynolds-number  flows. 
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Figure  7.1  Numerical  simulation  of  experimental  observations  by  Tollner, 
Barfield  and  Hayes  (1 982) 
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Vertical  Profiles  of  Suspended  Sediment 
Concentration 


After  the  overall  performance  of  the  numerical  code  was  checked,  vertical 
profiles  of  suspended  sediment  concentrations  were  computed.  It  was  found 
that,  due  to  the  particular  shape  of  the  vertical  eddy  viscosity  profile,  simulated 
relative  distributions  of  suspended  sediment  concentrations  differ  only  slightly 
from  the  classical  profile  in  open  channels  (i.e.  Rousean  distribution).  Figure  7.2 
and  7.3  illustrate  comparisons  between  simulated  dimensionless  profiles  of 
suspended  sediment  concentration  and  predictions  by  the  Rousean  model  for 
several  sediment  sizes,  and  two  different  conditions  of  bed  slope  and  plant 
density.  It  is  worth  noticing  that  relative  suspended  sediment  concentrations  are 
larger  than  the  ones  predicted  by  the  Rousean  formula.  It  would  seem  that  the 
effect  of  the  vegetation  is  to  promote  a  more  uniform  distribution  of  the 
suspended  sediment,  particularly  within  the  plants. 

Suspended  Sediment  Transport  Capacity 

Dimensional  analysis  has  shown  that  the  suspended  sediment  transport 
capacity  of  vegetated  channels  depends  upon  several  flow,  sediment  and  plant 
parameters.  This  information  is  herein  used  to  design  a  set  of  numerical 
experiments  aiming  at  characterizing  the  particular  effects  of  each  of  these 
parameters  on  the  transport  capacity.  Such  exercise  provides  also  a  means  of 
testing  the  reliability  of  the  numerical  results,  since  the  transport  capacity  should 
remain  the  same  when  dimensional  variables  are  changed  for  identical  values  of 
the  dimensionless  parameters.  Finally,  the  possibility  of  collapsing  the 
information  on  transport  capacity  for  different  conditions  is  also  investigated. 

Suspended  Sediment  Transport  Capacity  as  function  of  H/hp 

In  order  to  investigate  the  effects  of  the  ratio  H/hp  upon  the  suspended 
sediment  transport  capacity,  the  numerical  model  was  used  with  constant  values 
of  mean  flow  depth  and  plant  density  (H  =  0.35  m  and  a  =  2.0  m“0  while  three 
different  values  of  plant  height  were  used,  namely  hp  =  0.05,  0.10  and  0.25  m  in 
conjunction  with  varying  sediment  sizes  and  channel  slopes.  Figure  7.4  depicts 
the  results  obtained.  The  suspended  transport  capacity  has  been  scaled  using 
M*  Dj  (see  Yalin  1977);  however  any  other  scaling  may  be  easily  obtained 
combining  the  dimensionless  parameters.  For  example  Einstein’s  dimensionless 
version  becomes: 

^ss  _  9ss _ 

with  M*/  yjg  R  Ds  being  the  square-root  of  the  Shields’  stress. 
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k-e  Model 
Rousean  Model 


So  =  0.05  a  =  2.50m-^ 
H  =  0.35  m  hp  =  0.1175  m 


Dsf,  Dss  Ds4  Dsi  Ds2  Z)i, 


0.0  0.2  0.4  0.6  0.8  1.0 


c/c. 


Figure  7.2  Observed  and  predicted  vertical  distribution  of  dimensionless 
suspended  sediment  concentration,  condition  I 
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k-e  Model  So  =  0.0036  a  =  1.09 

Rousean  Model  H  =  0.35  m  hp  =  0.1175  m 


c/c. 


Figure  7.3  Observed  and  predicted  vertical  distribution  of  dimensionless 
suspended  sediment  concentration,  condition  II 
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Figure  7.4  Computed  suspended  transport  capacity  as  a  function  of  u-/Ws,  Rep 
and  H/hp,  for  H  a  =  0.70 
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It  was  observed  that,  as  expected,  transport  capacity  wtis  a  function  of  both 
m/ws  and  Rgp  for  each  value  of  the  ratio  H/hp,  so  that: 


Qss 

u*  Ds 


(94) 


where  the  function  /  will  depend  on  H/hp  (and  Ha).  Using  this  approach,  the 
data  could  then  be  collapsed  onto  a  single  curve  for  each  H/hp  ratio.  Results  are 
shown  in  Figure  7.5.  Moreover,  since  all  curves  in  the  previous  figure  are 
parallel,  it  is  easy  to  make  them  all  collapse  into  one  single  curve  writing: 


Qss 

u*  Ds 


=  f 


u* 

Ws 


(95) 


Figure  7.6  illustrates  results  for  the  best  fit  value  of  c  =  0.38. 


Suspended  Sediment  Transport  Capacity  as  function  of  Ha 

Following  similar  procedures  as  before,  values  of  mean  flow  depth  and  plant 
height  were  kept  constant  in  the  numerical  study  (H  =  0.25  m  and  hp  =  0.071  m), 
while  plant  density  was  allowed  to  vary,  namely  a  =  0.5,  2.8  and  5.0  m~^ 

Notice  that  for  a  =  2.8  m~i,  H/hp  =  3.5  and  a  =  0.7,  which  corresponds  to  one 
of  the  cases  studied  in  the  previous  section,  namely  in  Figure  7.4  b).  Since  in 
both  cases  values  of  these  two  parameters  were  equal,  although  the  value  of  each 
dimensional  variable  was  different,  similar  relations  were  expected  then  between 
the  dimensionless  transport  capacity  and  m/ws  for  each  Rep.  Results  are 
illustrated  in  Figure  7.7,  which  shows  a  good  collapse  of  the  computed  curves 
for  both  cases.  Figures  7.8  depicts  estimated  values  of  dimensionless  transport 
capacity  as  a  function  of  both  u*/ws  and  Rgp  for  each  value  of  the  parameter  H  a 
(H/hp  was  kept  constant  and  equal  to  3.5). 

Relative  Transport  Capacity  of  Suspended  Sediment 

The  influence  of  vegetation  in  reducing  the  suspended  transport  capacity  of  a 
channel  was  studied  by  computing  the  ratio  qs-veg^^s-oc  as  a  function  of  plant 
density  for  five  different  sediment  sizes.  Both  water  discharge,  q^,  and  plant 
height  were  kept  constant .  Here  qs-veg  is  the  computed  suspended  sediment 
transport  capacity  of  the  vegetated  waterway  and  qs-oc  is  the  capacity  for  an 
open  channel  of  the  same  slope  without  vegetation  and  same  water  discharge. 
Results  are  illustrated  in  Figure  7.9.  Note  that  in  keeping  q^  and  hp  constant 
both  H  a  and  H/hp  were  allow  to  vary.  It  is  interesting  to  note  that,  for  very  low 
densities,  the  ratio  qs-veg^qs-oc  becomes  slightly  larger  than  unity  for  the  coarsest 
sediment.  This  apparent  contradiction  may  be  explained  by  the  fact  that  (at  low 
densities)  a  small  increase  in  density  tends  to  decrease  the  momentum  transfer 
toward  the  bed  only  by  a  small  amount  while  at  the  same  time  requiring  a  small 
increase  in  flow  depth  due  to  the  increased  flow  resistance.  The  combined  effect 
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Figure  7.5  Computed  suspended  transport  capacity  as  a  function  of  the 
parameters  u»Ms  flgp^and  H/hp,  for  Ha  =  0.70 
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Figure  7.6  Computed  suspended  transport  capacity  as  a  function  of  parameter 
U-/Ws  Rep‘='  (H/hp)°-^^  for  Ha  =  0.70 
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Figure  7.7  Computed  suspended  transport  capacity  as  a  function  of  U’/Wg  and 
Rep’  for  H  a  =  0.70  and  H/hp  =  3.5.  Closed  symbols  correspond  to 
{H,a,hp)  =  (0.35  m;  2.0  m"'':  0.10  m),  whereas  open  symbols  are  for 
(H,a,hp)  =  (0.25  m;  2.8  m-i;  0.07  m) 
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Figure  7.8  Computed  suspended  transport  capacity  as  a  function  of  u^/Wg,  Rgp 
and  Ha,  for  H/hp  =  3.5 
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tends  to  keep  the  shear  stress  at  the  bottom  nearly  constant,  and  therefore  to 
entrain  as  much  sediment  from  the  bed  as  without  the  plants.  Moreover,  if  it  is 
remembered  that  the  relative  sediment  concentration  profile  shows  larger  values 
in  the  presence  of  vegetation,  it  can  easily  be  explained  why  suspended  transport 
capacity  may  become  slightly  larger  in  the  presence  of  plants  compared  to  flow 
without  vegetation. 

Final  Remarks 

Dimensional  analysis  helped  in  identifying  the  different  dimensionless 
parameters  that  govern  sediment  transport  processes  in  vegetated  water  channels, 
and  the  numerical  model  proved  to  consistently  predict  suspended  sediment 
loads  under  different  conditions  for  the  same  parameter  values.  It  is  worth 
mentioning  that,  albeit  at  a  preliminary  level,  a  different  type  of  two-equation 
model  has  also  been  developed  as  an  alternative  code,  namely  a  k-w  type  closure 
(Lopez  and  Garcia  1996).  Results  of  both  models  have  been  found  to  provide 
similar  degree  of  representation  to  the  experimental  observations.  Only  as  an 
example.  Figure  7.10  shows  the  dimensionless  sediment  transport  capacity  as 
computed  by  the  k-oo  model,  compared  to  the  best  fit  of  results  obtained  with  the 
k-E  model. 
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Figure  7.10  Suspended  sediment  transport  capacity  as  a  function  of  parameter 
u»/Ws  Rep^  (H/hp)^-^^,  for  /-/a  =  0.70.  Computations  by  the  k-co  model 
compared  to  best-fit  line  to  results  from  k-e  model. 
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8  Summary  and  Conclusions 


From  the  results  shown  in  the  previous  sections,  the  following  conclusions 
can  be  drawn: 

a.  The  decrease  in  suspended  sediment  transport  capacity  of  channels  with 

vegetation  compared  to  non- vegetated  ones  at  similar  water  flow  rates  (i.e. 
the  retention  capability  of  the  former)  is  highly  dependent  on  the  reduced 
ability  of  flow  over  a  vegetation-covered  bed  to  entrain  sediment  in 
suspension.  This  fact  is  in  turn  associated  with  the  decrease  in  streamwise 
momentum  transferred  to  the  channel  bed,  due  to  the  absorption  of 
momentum  by  the  plants  via  drag  forces. 

b.  The  above-mentioned  fact  has  two  important  practical  consequences,  the 
first  one  being  that  any  simpler  model  developed  to  compute  the  transport 
capacity  of  vegetated  channels  should  be  based  on  a  reliable  estimate  of 
the  average  shear  stress  taken  by  the  bed,  which  implies  a  good 
characterization  of  the  form  drag  coefficient  of  plants  in  water  channels. 
The  second  one  is  that  laboratory  and  field  studies  are  needed  in  order  to 
develop  a  sediment  entrainment  function  for  flow  through  vegetation.  The 
characteristics  of  the  near-bed  turbulence  in  vegetated  open-channel  flows 
are  such  that  wake-generated  turbulence  might  be  the  main  mechanism 
responsible  for  sediment  entrainment  into  suspension,  playing  a  role 
similar  to  that  of  turbulent  bursts  in  boundary-layer  flows  without 
obstructions. 

c.  The  two-equation  model  of  turbulence  developed  herein  provides  a  good 

representation  of  the  experimental  observations  of  different  turbulence 
variables,  length  scales  and  energy  budget  terms,  hence  constituting  an 
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alternative  tool  for  analyzing  the  influence  of  different  flow  and  vegetation 
properties  on  the  overall  turbulence  structure. 

d.  The  numerical  model  consistently  predicts  the  sediment  transport  capacity 
under  different  flow,  sediment  and  vegetation  conditions  for  same  values 
of  the  governing  dimensionless  parameters,  and  proper  combination  of 
these  parameters  further  allows  for  the  collapse  of  all  information  onto  one 
single  relation. 

e.  The  Rousean  profile  of  relative  sediment  concentrations,  computed  with  a 

shear  velocity  derived  from  the  total  action  of  gravity  forces  (i.e. 

M*  =  -Jg  H  So),  predicts  relative  distributions  very  similar  to  the  ones 
obtained  with  the  numerical  model.  This  might  be  mainly  attributed  to  the 
parabolic-type  shape  of  the  eddy  viscosity  profile  for  both  open  channels 
with  and  without  vegetation. 

f.  Model  results  show  the  Manning’s  coefficient  to  remain  almost  constant 

(with  values  close  to  non-vegetated  conditions)  up  to  a  critical  plant 
density,  and  to  increase  linearly  afterwards,  in  agreement  with  field 
studies. 

g.  In  summary,  the  two-equation  numerical  code,  based  on  the  k-e  turbulence 
closure  scheme,  constitutes  a  reliable  tool  for  engineering  assessments 
both  on  the  turbulence  stmcture  and  the  suspended  sediment  transport 
processes  in  open  channels  through  vegetation.  The  challenge  for  the 
future  consists  in  extending  the  capabilities  of  the  model  developed  for 
“idealized”  vegetation  to  the  case  of  natural  plants. 

h.  Future  numerical,  laboratory,  and  field  work,  will  also  benefit  from  the 
dimensional  analysis  carried  out  in  the  previous  chapter.  Also  of 
particular  relevance  for  future  laboratory  and  field  measurements,  is  the 
analysis  conducted  in  Chapter  3,  where  the  need  to  perform  spatial  as  well 
as  temporal  averaging  of  flow  measurements  in  order  to  obtain  meaningful 
results  was  clearly  demonstrated. 


Chapter  8  Summary  and  Conclusions 
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Appendix  A  Notation 


(.)  Time-averaged  operator  over  turbulence 

<•>  Spatially  averaged  operator 

(.)'  Fluctuation  over  time-averaged  value 

(.)"  Fluctuation  over  space-averaged  value 

(.)+  Variable  made  dimensionless  using  wall  units  (U*  and  v) 

V  Laplacian  operator 

a  Ratio  between  the  sum  of  the  differential  frontal 
areas  of  the  obstacles  divided  by  the  differential 
volume  of  fluid 

a  Dimensionless  parameter  defining  flexibility  of  vegetation 

b  Distance  equal  to  five  percent  of  the  flow  depth  measured 

from  the  bottom 
bp  Plant  spacing 

B  Bouyancy  parameter  5  =  -  gR  (  dC/dz  ) 

C  Suspended  sediment  concentration 

Q?  Bottom  sediment  concentration 

Cc'  Model  parameter  for  temperature  stratification  effects  on 
Prandtl-Schmidt  number 
Q)  Drag  coefficient  of  plants 

Cds  Drag  coefficient  of  sediment  particles 

Eddy- viscosity  coefficient  in  the  k-e  turbulence  model 
Co  Averaged  suspended  load  concentration 

C],  C2  Weighting  coefficients  for  production-  and  dissipation-related 

terms 

C]R  Pressure-strain  coefficient 

Q  Coefficient  modifying  the  traditional  wall  function  for  s 

Ck  Coefficient  modifying  the  traditional  wall  function  for  k 

Cvei  Coefficient  modifying  the  traditional  wall  function  for  U 


Appendix  A  Notation 


A1 


Cfs 

Weighting  coefficients  for  production-  and 
dissipation-related  drag-terms 

D 

Horizontal  diameter  of  the  (vertically  oriented)  plants 

Ds 

Mean  sediment  diameter 

^ij 

Kronecker  delta 

At,  Az 

Time  and  vertical  spacing  on  the  numerical  grid 

E 

Roughness  parameter  approximately  equal  to  9 
for  hydraulically  smooth  conditions  and  30v/u*/ks 
for  fully  rough  conditions  where  kg  =  equivalent 
sand  roughness 

El 

Stem  flexural  rigidity 

e 

Dissipation  rate  of  turbulent  kinetic  energy 

Fr 

Froude  number 

fx 

Pressure  Force  per  unit  length  in  z  on  the  perimeter  s 

fi 

Drag  force  per  unit  volume,  i-direction  component 

0c>  ^c' 

Model  parameters  for  temperature  stratification  effects  on 

Prandtl-Schmidt  number 

<Ps^ 

Modified  version  of  Grass’  parameter 

8 

Gravitational  acceleration 

8i 

Gravitational  acceleration  component  in  the  i-direction. 

Y 

Pressure-strain  coefficient 

hp 

Plant  height 

H 

Mean  flow  depth 

V 

Kolmogorov  microscale 

kc 

Measure  of  roughness  in  the  Einstein’s  (1950) 
equation  for  suspended  load 

ke 

Einstein’s  viscosity  constant 

ks 

Equivalent  sandgrain  roughness 

K 

Von  Karman’s  constant  (0.40) 

kw 

One-dimensional  wave  number  in  the  streamwise  direction 

k 

Turbulent  kinetic  energy 

ko 

Turbulent  kinetic  energy  evaluated  at  the  first  grid  point 
from  the  bed 

1 

Mixing  length 

h 

Length  scale  for  momentum  transfer 

4 

Macro-length  scale  for  streamwise  velocity  fluctuations 

X 

Taylor  microscale 

4 

Plant  density  in  stems  per  square  meter 

M 

Relative  density  of  plants 

A2 

Appendix  A  Notation 

N  Total  number  of  grid  points  in  the  vertical  direction 
n  Manning’s  resistance  coefficient 

”  Vector  normal  to  the  perimeter  of  the  object 
Component  of  n  in  the  x-direction 

>^xu’  f^xd  Vector  nx  in  the  upstream  and  downstream  faces  of 

the  object,  respectively 
V  Fluid  kinematic  viscosity. 

^  Kinematic  viscosity  of  sediment-water  mixture 
Kinematic  viscosity  of  clear  water 
vr  Kinematic  eddy  viscosity 

Fluid  dynamic  viscosity 

P  Instantaneous  pressure 

p-^  Shear-  and  wake-production  terms 

W  Modified  version  of  Grass’  parameter 

Q  Total  water  discharge 

qw  Specific  water  discharge 

qss  Suspended  sediment  transport  capacity 

qs-oc  Suspended  sediment  transport  capacity  without  vegetation 

qs-veg  Suspended  sediment  transport  capacity  with  vegetation 

Re  Flow  Reynolds  number 

Rep  Dimensionless  particle  size  defined  as 

Rep  ~  '] 8  R  f 

Rf^  Hydraulic  radius 

Rs  Equivalent  hydraulic  radius 

R  Submerged  specific  gravity  of  sediment,  defined  as 

=  (Qs-Q)/q 
Q  Fluid  density 

Qs  Sediment  density 

Pw  Clear  water  density. 

s  Perimeter  of  a  cylinder 

Generalized  source  term 
Bed  slope 

S(kw)  One-dimensional  normalized  spectra  for  streamwise 
velocity 

a  Prandtl-Schmidt  number 

Prandtl-Schmidt  number  for  sediment  particles 
Oco  Prandtl-Schmidt  number  for  non-stratified  conditions 

t  Time 
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Tt 

Ts 

U]  (u),  U2(v),  U3(w) 

U,V,W 

/  /  / 
U  ,  V  ,  w 

^rms*  ^rms 

Uo 

m 

^*hp 

x,y,  z 

Ws 

Zu 

Zo 

X 


Turbulent  transport  of  turbulent  kinetic  energy 
Time-scale  of  momentum  transfer 

Instantaneous  streamwise,  spanwise  and 
wall-normal  velocities,  respectively 
Mean  streamwise,  spanwise  and  wall-normal 
velocities,  respectively 

Streamwise,  spanwise  and  wall-normal  velocity 
fluctuations,  respectively 

Root-mean  square  values  of  streamwsie,  spanwise  and 
wall-normal  velocity  fluctuations,  respectively 

Value  of  U  at  first  grid  point  from  the  bed 
Mean  bed  shear  velocity 

Square-root  of  the  Reynolds  stress  per  unit  density 
at  the  top  of  the  canopy 
Right-handed  coordinate  system  representing 
streamwise,  spanwise  and  wall-normal  axis,  respectively. 

Terminal  fall  velocity  of  sediment  particle 

Garcia  and  Parker’s  parameter  {Zu  =  m*/^^  R^p  ) 

First  grid  point  away  from  the  bed 

Generalized  variable  characterizing  sediment  transport 

process 
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